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SUMMARY 


The  present  review  contains  a  description  with  comments  of  the 
three  main  Monte  Carlo  methods  which  have  been  used  to  date  in  solution 
of  the  full  Boltzman  equation  of  kinetic  theory  -  namely 

lf  Test  particle  method  of  J.  K.  Haviland. 

2.  Simulation  method  of  G.  A.  Bird. 

3.  Integral  Evaluation  method  of  A.  Nordsieck  and  B.  L. 

Hicks . 

A  chapter  containing  the  necessary  formulas  from  kinetic  theory 
and  one  on  probability  theory. is  included  at  the  beginning.  The  author's 
first  hand  experience  with  the  simulation  method  has  made  possible  the 
inclusion  of  some  new  material  in  Chapter  U. 
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1.  PRELIMINARY- NOTATION-MD'THEORY 


1 • 1  Elastic  Collision  of  Two 'Molecules  Having  Equal  Mass. 

1.1.1  Classical . Description 

A  Collision  between  two  such  molecules  M  and  M'  is  adequately- 
described-  in  this  setting  whenever -two  pre-collision  velocity  vectors 
v,  v'- and"two  postrcollision  velocity  vectors  V,  V'  are  giyen. -  The 
terms  pre-collision  and  post-collision  are  only  meaningful  when  the 
molecules  move  freely  before  and  after  collision,  that  is  when  the  forces 
between  them  act  only  when  they  are  close  together.  This  assumption  is 
necessary  to  even  speak  of  a  collision  process. 

The  theory  shows  that  V,  V"  are  fuuctiqns  of.  v,  v'  and  two  other 
parameters  b,  e  whose  geometrical  meaning  is  explained  in. Fig.  1  (b  is 
called  the  miss  distance  and  e  specifies  the  orientation  of  the  collision 
plane).  Figures  1-3  refer  to  a  coordinate  system  moving  with  velocity  v 
in  which  M  is  initially  at  rest  at  the  origin.  The  pre-collision  and 
postrcollision  relative  velocity  vectors  are  denoted  by  vrj  =  v'-v  and 
Vr  =  V'-V  respectively  and  x  is  the  deflection  angle  of  M.  It  can  be 
shown  that  Vr  lies  in  the  collision  plane  (vr  lies  in  this  plane  by 
definition,  so  Vr,  vr  are  coplanar),  and  that  x  =  x(b,g)  (g  Hvr|) 
where  the  function  x  is  determined  by  the  collision  interaction.  For 
purposes  of  kinetic  theory  it  is  convenient  tp  assume  that  there  is  a' 
value  ,bm  such  that  x  =  0  for  b  >  bm  as  would  be  the  case  for  a  purely 
local  interaction.  It  is  still  a  very  good  approximation  in  cases  of 
physical  interest. 


The  collision  vector  e  in  Fig.  3  is  a  unit  vector  in  the  direction 
of  momentum  transfer  during  collision  and  is  completely  determined  in  the 
collision  plane  by  b  and  Vr,  and  in  three  dimensions  by  v,  vr  and  e 
(e  =  e(vr,b,e).  Figure  k  is  the  analogue  of  Fig.  3  for  an  arbitrary 
coordinate  system  in  which  M  and  M'  have  pre-collision  velocity  vectors 
v  and  v"  respectively.  The  following  formulas  for  the  functions  V  and 
V"  are  useful  for  subsequent  work: 


V  =  v  +  (vr  •  e)e 
V"  =  v”—  (vr  •  e)e 


(1) 


V  =  l/2(v"  +  v  -  Vr) 

(2) 

V-  =  l/2(v'  +  v  Vr) 

1.1.2  Statistical  Description  -  Scattering  Cross-Section 

The  basic  element  in  the  statistical  description  of  the  collision 
process  is  the  conditional  probability  density  (differential  cross-section) 

G (tf,  g)  having  the  property  that  G($,  g)  d#  is  the  probability  that  the  post 
collision  velocity  vector  of  M"  has  direction  within  a  solid  angle  dO  about 
the  unit  vector  $  given  that  a  collision  has  takon  place. 
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G  can  "be  obtained  in  terms  of  X  by  first  observing  that  the 
distribution  for  (b,e)  given  b  <_  bm  has  the  density  function  b/ub^. 

The  calculation  is  simple  if  we  assume  that  an  inverse  function 
h(X,g)  can  be  found  with  the  property  that 

X(b(X,g),g)  =  X 

Then  the  probability  that  X  lies  in  the  range  (x,X+dX)  and  e  lies  in  the 
range  (e,e+de)  is  given  by 


TT*  !£(x’g)  dXde 

m 

Since  the  range  (X,X+dX)  x  (e,e+de)  corresponds  to  a  solid  angle  sin  x^xde, 
it  follows  that 


where 


G(w,g)  = 


3b 

irb  ^siny  3y 
m 


(x,g) 


w=  (cosxcose,  cosxsine,  sinx) 


(3) 


1 . 1 i 3  '  The  Hard  Sphere  Model 

This  model  constitutes  one  of  the  rare  cases  for  which  x  and  G 
have  explicit  expressions  and  is  the  one  most  often  used  in  Monte  Carlo 
calculations.  The  key  formula  for  our  purposes,  evident  from  Fig,  5  is 

b  cos  x  =  b. 
m  A 


It  follows. that 


db  =  -  b^  sinxdx, 

and  the  probability  that  (x,e)  lies  in  (x»X+dx)  x  (e,e+de)  is 


b^  cosx'  bm  sinx  dx  ds 


nb  ^ 


m 


sin  2x  d(2x)  d£ 
Utt 


Referring  again  to  Fig.  5,  we  see  that  if  the  deflection  angle  is  x  then 
the  change  in  direction  of  the  relative  velocity  vector  is  2x.  Now  the 
solid  angle  corresponding  to  the  range  (2x,  2x  +  d(2x)  (e,e  +  de)  is 

sin  2x  3i('j2x)  de. 


20. 


FIGURE  5:  HARD  SPHERE  MODEL 


From  the  preceeding  two  facts  we  deduce  that  Vr  is  uniformly  distributed 
on  the  surface  of  the  unit  sphere,  and  most  surprisingly,  distributed 
independently  of  vr.  This  fact  makes  the  calculation  of  post-collision 
velocities,  very  simple  in  the  hard  sphere  model. 


1.2  The  Boltzman  Equation 

•  1,2,1  The  Distribution  Function  vf 

Suppose  we  wish  to  give  a  non-aetailed  description  at  time  t 
of  a  gas  of  like  molecules  in  some  region  <0  of  three  dimension  Euclidean 
space  E3.  One  way  of  doing  this  is  to  specify  a  positive  function 
f(x,V,t)  for  ice  £>  ,  veE3  such  that  f(x,v,t)  dx  dv  yields  the  expected 
number  of  molecules  in  the  small  six-dimensional  rectangular  cell  about 
the  point  (x,v)  whose  sides  have  lengths  given  by  (dx,dv).  The  set  of 
all  pairs  (x,v)  is  known  as  phase  space  and  is  denoted  by  fi  =  <D  x  E3, 

1.2.2  The  Equation  for  the  Classical  Description 

The  Boltzman  equation  is  a  conservation  equation  for  f.  Its 
basic  form  is 


(||-  +  v.Vf)  (x,v,t) 

=  J  dv'  bdbJo  delvJf^V't )f(£,V,t)  (h) 


-J  dv'J  m  bab  de jvr|f(x,v',t)f(x,v,t) 

where  V'  and  V  are  functions  of  v',  v,  b,  s  as  Ascribed  in  the  foregoing 
section.  The  three  terms  which  have  been  displayed  in  the  above  equation 
are  known  respectively  as  the  convection  term,  the  gain  term,  and  the  loss 
term. 

Because  the  expressions  involved  are  lengthy  ones,  certain 
conventions  have  been  adopted  in  order  to  snorter)  the  amount  of  writing 
involved,  Let  f  denote  f(x,v,t),  f'  denote  f(x,v'.t),  F  denote  f(x,V,t), 
and  F'  denote  f(x,V',t),  Then  (k)  can  be  written  as 

a-r  r  z'*’  r2n 

~  +  v  •  Vf  =  j  dv'  j  m  bdb  j  de|v  jv^F-f'f)  (5) 

An  examination  of  the  loss  term  reveals  that  f(x,v,t)  is  independent  of 
the  variables  of  integration  and  can  therefore  be  taken  outside  the 
integral  sign.  Accordingly  the  gain  term  is  often  denoted  by  A  and 
the  loss  term  by  Bf  .  Thus 
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(6) 


77-  +  v  *  w*  A  -  Bf 

(it 


The  combination  A-Bf  is  often  referred  to  as  the  collision  integral  and 
written  C(f)  to  emphasize  that  it  is  an  operator  on  f. 

The  mathematical  representation  of  a  physical  problem  requires 
that,  in  addition  to  equation  (4)  initial  and  boundary  conditions  on  the 
function  f  must  be  satisfied.  If  the  .representation  is  complete  it  is 
then- expected  that  the  solution  f  will  exist  and  be  unique . 


1.2.3  Simplification  for  Hard  Spheres 

If  we  recall  that  b  =  bjflCosx  (see  Fig.  5)  where  x  is  the  angle 
between  vr  and  if,  it  follows  that 


|vr|b  db  de  =  b^2  (if  •  vr)  dfr 
Consequently  (?)  becomes 

||  +  v  •  Vf  =  b m2f  dv'  J  d£_  (if  •  vr)  (F'F-f'f) 

°  "  #  *  Vr  >_  ° 

and  by  making  a  minor  extension  in  the  definition  of  V  and  V'  this  can  be 
written  as 


3f 

3t 


+  v 


Vf  =  1/2  bjj yj  dv'  f  d^l'fr  •  vr|  (F'F-f'f) 


(7) 


1.2.h  The  Equation  in  Terms  of  G(fl,g) 

The  Boltzman  equation  can  also  be  written  in  terms  of  the 
Ifferential  cross-section  G(if, g).  The  expressions  in  this  case  are 


(||  +  v  •  Vf)  (x,  v,  t) 


(8) 


=  J dv'  J  dtfjvr|G(tf.|vr|)  (f(x,  V',  t)f(x,  V,  t)-f(x,  v',  t)f(x,v,t)) 


where  now,  of  course,  V,  V"  must  be  given  as  functions  of  v,  v' ,  if,  instead 
of  v,  v',  b,  e.  The  integral  over  if  in  this  equation  is  non-zero  only  over 
the  forward  hemisphere  0  <  x  <  n/2  since  G  =  0  for  u/2  <  x  <  it  to  allow 
conservation  of  momentum. 


It 


1.2.5  Cylindrical  Coordinates  in  Velocity  Space 


In  problems  requiring  only  one  space  dimension  to  describe,  The 
distribution  function  f  often  has  tee  form 


f  (x,y-z,v;,vyavz,t)  =  ffx,’^,^2^2)27* 
=  f(x,T ■  &_  ) 


(9) 


where  v,  =  (v  2  +  v  2)l/2 

y  2 

In  this  case  the  left  side  of  (U)  has  the  form 


of  A  3f 
3t  Vx  3x 

Although  no  reduction  in  the  dimension  of  the  collision  integral  follows, 
it  is  fairly  easy  to  show  that  for  any  f  of  this  fora,  C(f)  has  the  seme 
form. 


To  this  end  it  is  useful  to  introduce  the  representations 


v  =  (v  ,  t^cos^,  sin p),  V  =  (v/,  vx  "  cosf  %  vx  'sin p  '), 

and 

w  =  (wx,wxcosV,  v^sinv)§  (10) 

where 

(wx,wx)  =  (cos  0,  sin  0), 

and  to  write  the  collision  integral  with  respect  to  dvx  dv^  djp'  d0  dv. 

It  turns  out  that  the  integrand  in  C(f)  depends  only  on  the  differences 
[lp-v)  and  and  not  on  the  individual  angles.  Since  the  integral? 

with  respect  to  ip'  and  dv  are  over  the  interval  (0,  2ti),  a  simple  changr 
of  variables  yields  the  result  that  C(f)  is  independent  ofjp  .  The  Boltzc  jn 
equation  using  cylindrical  coordinates  in  velocity  space  is  then 

3f  3f  \2  ra  r"  f2tr 

St  *vxl t’-T  }  aVJ  'i'  aV J  af  a5|w  •  ; 

J  — cc  o  O 

(n) 

where  the ’’'s  have  been  omitted  for  convenience. 


1.3  Invariance  Properties 

1.3.1  Normalization  of  f 

^  If  f  satisfies  equation  (8)  with  differential  cross-section 

G(w,g)  and  a  >  0  then  af  satisfies  (5)  with  differential  cross-section 
G(w?g)/a,  This  makes  it  possible  to  work  with  a  normalized  distribution 
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function  f  having  the  property  that 


v,t  )  av  =  1 
o 


at  some  reference  noint  x  and  tine  t  .  That  is. 

o  o 


' 

=  f(x,vst)/no 

where 

r 

n  =  j 

f(x  ,v,t  )  dv 

°  JE3 

o  o 

The  new  differential  cross-section  is  then 


C-(tf,g)  =  nQ  G(*,g) 


and  can  he  interpreted  as  a  differential  cross-section  density  because  of 
the  equation 


C-(tt,g)  dtt 


-O  S'**’ 


g)  d$- 


n  A(g). 
o 


The  quantity  (42  n  A(g))'1  is  the  nean  free  path  at  speed  g.  In  the 
case  of  hard  sphere  collisions  both  A  and  the  mean  free  path  are  indepen¬ 
dent  of  g.  In  this  case,  if  L  is  a  characteristic  dimension  in  the 
problem. 


L  (nA)-1 

constitutes  a  useful  dimensionless  matching  parameter  commorjly  known  as 
the  Knudsen  number. 


1.3.2  Change  of  Scale 

In  theoretical  work  it  is  frequently  convenient  to  express 
the  basic  equations  in  terms  of  time  and  length  units  characteristic 
of  the  problem  at  hand.  To  this  end  let 

t  =  Tt1} 
x  =  A.^, 
and  a  =  Vx, 

from  which  it  follows  that 

V  =  ovr 
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Put 


Then 


f(xj  ,ti)  =  (Ac)3f(Axi,  jCrvjjTtj). 


f  f(r2,v i,t5)  dx2  dv2  =  (  f(x,v,t)  dx 

Furthermore,  if  f (x,v,t)  satisfies  equation  (k)  with 


av. 


V  =  v  -r  (v  -e)  e 


V'  =  v'-(v-e)  e 


e  =  e  (v  ,hse)s 


then  a  straightforward  calculation  reveals  that  f(x,v,t)  satisfies. 

-2w 


3f 

3t2 


b  /X 

•  El 


vi-Vif  =  fe'f  b2db2J  ds | Vjr | (?'F-£'t)  (12) 


vnere 


f  =  fCxjjVjjtj) 

f'= 

F  =  f(x,  ,?!*!) 

F  =  ftx^VjStj) 
V2= 

V^l  "+(vlr-e1)e1 


That  is,  the  Maxwell-Bolt zman  equation  is  invariant  imder  scale  changes 
if  the  proper  change  is  made  in  the  collision  vector  5  and  the  miss 
distance  h^. 

l.U  The  Boltzaan  Equation  from  a  Computational  Standpoint 

Any  attempt  to  solve  the  full  Boltzman  equation  by  computational 
means  must  offer  a  solution  to  two  basic  problems  -  a  feasible  evaluation 
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of  the  collision  integral  and  a  convergent  strategy  for  obtaining  a 
solution  to  the  first  order  nonlinear  partial  differential  integral 
equation  which  constitutes  the  Boltzman  equation.  A  straightforward 
resolution  of  these  problems  fails  because  both  the  fast  access  storage 
capacity  and  the  computing  speed  of  even  our  fastest  machines  are 
inadequate  for  dealing  with  even  simple  problems.  There  are,  of  course, 
other  questions  to  be  considered  such  as  boundary  conditions,  but  they 
pose  difficulties  of  a  lesser  order  of  magnitude. 

Even  the  storage  requirements  for  the  distribution  function 
f(x,v,t)  at  fixed  time  t  are  considerable,  for  suppose  that  f  were 
tabulated  by  storing  its  values  at  10  grid  points  for  each  degree  of 
freedom.  Then  in  the  general  case  of  3-dimensional  x  and  v,  106  fast 
access  storage  locations  would  be  required  which  is  beyond  the  capacity 
of  present  day  computers.  In  a  highly  symmetrical  problem  such  as  the 
normal  plane  shock,  one  can  get  by  using  one  dimension  for  x  and  two 
dimensions  for  v  by  making  use  of  the  symmetry.  In  this  example  only 
103  values  are  required  under  the  assumption  of  10  grid  points  per 
degree  of  freedom.  Since  103  is  well  within  the  range  of  most  computers, 
this  particular  problem  has  been  done  by  several  authors. 

The  evaluation  of  the  collision  integral  by  numerical  quadrature 
is  considerably  more  intractable.  To  gain  an  appreciation  of  the  diff¬ 
iculties  involved  it  is  sufficient  to  consider  a  typical  calculation. 

Since  x  and  t  remain  fixed  during  one  evaluation  of  the  collision  integral, 
the  distribution  function  will  be  written  f(v);  and  since  it  is  expected 
that  f  decreases  very  rapidly  (like  exp(-c|v|2))  for  large  v,  consider 
f  only  in  a  subregion  S  of  velocity  space  with  the  property  that  the 
integral  of  f  outside  of  S  is  negligible.  Suppose,  furthermore  that  S 
is  subdivided  into  bins  denoted  by  S^  1  <  k  <  n,  and  that  f  will  be 
represented  in  computer  memory  by  storing  the  n  numbers  f^= /g.  f  (v)dv 
1  k  <_  n.  That  is,  the  actual  f  used  in  the  computation  will  be  an 
approximation  given  by 


-  L  f  A(">  <13) 


where  y.  (v)  is  a  function  taking  the  value  1  for  v  e  S,  and  0  otherwise. 
Then  k 

b  2tt 

C(f)  =  E  f.fj  Jdv'  J  bdb.J  de  [xi(V")xj(V)-X-(vOxj(v)]  (lU) 


and  the  resultant  contribution  to  bin  k  is  given  by 


°*(f)  =  ?  A 

1»J 


id 


f.f. 
1  J 


(15) 
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where 


b  2  77 

f  dv  fdv"  f  hdb  f  de  tx* (V')x*(V)-x.(v')x,(v)] 

JS  J  J  o  o  3  3 

fc 


Even  under  the  course  assumption  that  there  he  10  bins  in  each  of  three 
directions  in  velocity  space,  n  has  to  he  103.  This  means  that  the 
evaluation  of  C  requires  2  x  106  multiplications  and  106  additions 
at  each  point  k  or  2  x  109  multiplications  and  109  additions  altogether. 
Since  even  the  fastest  present-day  computers  have  multiplication  times 
of  the  order  of  10“6  sec.,  it  is  evident  that  capabilities  fall  short 
of  requirements  even  without  considering  the  colossal  task  of  computing 
and/or  storing  the  coefficient  array  P  . 

K 

The  problem  of  finding  the  best  strategy  for  solving  the 
nonlinear  differential  equation 


~  +  v  •  V  f  =  C(f ) 

is  at  least  partially  a  theoretical  one  and  presumably ^solvable  even 
in  the  absence  of  a  satisfactory  means  for  evaluating  the  collision 
integral.  Regretably,  very  little  theory  is  known  concerning  equations 
of  this  type.  At  least  two  methods  have  been  tried  in  the  case  when  C 
is  replaced  by  a  sampler  approximation  C^.  _The  first  is  concisely 
described  by  the  following  formula  for  i(x,v,t+At)  in  terms  of  f(x,v,t): 

f(x,v,t+At)  =  f(x,v,t)  +  [-v  •  Vf  +  C  (f)]  At.  (l6) 

Si  i 

In  the  case  of  time  dependent  problems  (l6)  yields  a  sequence  of 
approximations  to  the  distribution  function  at  intervals  of  time  At. 

In  the  case  of  time  independent  problems,  it  is  expected  that  the 
sequence  of  approximations  will  converge  to  a  solution  of  the  equation 

v  •  Vf  =  C  (f)  (IT) 

a 

in  much  the  same  way  as  a  non-equilibrium  distribution  function  relaxes 
to  its  equilibrium  steady  state  value  in  a  physical  situation.  This 
has  been  found  to  be  a  fairly  slow  method  and  greater  speed  is  achieved 
by  solving  (17)  directly  by  means  of  a  finite  difference  scheme. 

In  spite  of  the  major  hurdles  that  must  be  overcome,  a  tremen¬ 
dous  amount  of  work  has  been  dene  adding  to  our  knowledge  of  the  Boltzman 
equation.  A  whole  volume  would  be  far  from  sufficient  to  describe  it. 

The  major  categories  of  solution  methods  include  the  moment  methods,  use 
of  the  BGK  collision  integral,  and  most  recently,  the  Monte  Carlo  methods. 

In  this  review  we  describe  the  Monte  Carlo  methods  cf  J.  K. 
Havilland,  G.  A.  Bird,  and  A.  Nordsiek.  An  attempt  has  been  made  to 
include  a  large  number  of  related  works  in  the  bibliography. 

In  order  to  make  the  discussion  as  self  contained  as  possible, 
a  chapter  on  the  probability  theory  relevant  to  Monte  Carlo  calculations 
is  included  at  the  beginning. 
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2.  BASIC  STATISTICS  AND  PROBABILITY  THEORY 


2.0  Introduction 


The  purpose  of  this  chapter  is  to  assemble  some  of  the  ideas 
and  concepts  from  probability  theory  which  the  author  has  found  useful 
in  dealing  with  Monte  Carlo  methods.  No  attempt  shall  be  made  to  lay 
down  a  precise  enough  structure  to  permit  rigorous  proofs.  The 
definitions  are  more  complete .than  the  remainder  of  the  discussion  to 
permit  the  reader  to  easily  find  additional  information  in  a  probability 
theory  textbook,  but  in  most  cases  we  do  not  take  the  space  to  fully 
explain  all  the  hypothesis  mentioned  in  a  definition  if  they  are  commonly 
satisfied  in  the  usual  examples  of  probability-  spaces .  The  notation  in 
this  section  is  characteristic  of  probability  theory  and  independent  of 
that  used  up  to  this  point. 


2.1  Probability  Space 

2.1.1  Definition 


The  central  notion  in  statistics  and  probability  theory  is 
that  of  a  sample  space  or  probability  space,  and  we  begin  by  describing 
this  entity  in  very  general  abstract  terms.  Real  understanding,  of  course, 
comes  from  seeing  how  the  concept  is  used  in  a  number  of  familiar  situat¬ 
ions  and  we  shall  get  to  that  immediately  afterwards.  We  assume  for  the 
purposes  of  the  next  paragraph  that  the  reader  is  familiar  with  the 
concepts  of  set,  element,  subset,  union  and  intersection. 

Let  ft  be  a  set  and  V  be  a  o-field  of  subsets  of  The 

subsets  in  ¥  will  turn  out  to  represent  events.  The  term  o-field 
requires  that  ~f  contain  the  empty  set  <j>  and  that  it  be  closed  under 
the  operations  of  countable  union  and  set  difference.  Suppose,  in 
addition,  that  we  are  given  a  function  P  whose  domain  is  ¥  and  whose 
range  is  the  interval  [0,l]  satisfying  P(^» )  =  0,  P(ft)  =  1  and  p(v)  A  ) 

=  Z  P(A  )  whenever  the  A  ,  n  =  1,  2,....  are  disjoint.  This  n  ” 
last  condition  implies  tlmt  P(AUB)  =  P(A)  +  P(B)  -  P(AftB)  as  well  as  a 
certain  continuity  property  for  P.  We  call  the  triple  (ft,^  ,P)  a 
probability  space. 

2.1.2  Coin  Tossing 

We  will  try  to  show  that  the  abstract  structure  of  the  last 
paragraph  provides  us  with  a  model  fitting  all  sorts  of  situations 
which  according  to  our  intuition  are  of  a  probabilistic  nature.  Take, 
for  example  the  common  coin  tossing  experiment. 

Consider  first  a  single  toss  of  a  fair  coin.  We  claim  that 
(ft, ¥  ,P)  with  ft  =  Jh,tI,  where  H  denotes  heads  and  T  denotes  tails, 

*  =  ft,  {h}  ,  (t}\  {h,  T}}  and  P($)  =  0,  P(  {h}  )  =  P(  (T }  )  =  1/2, 

P(  {K,t}  J=1  mode],  for  this  experiment.  The  reason  it  is  suitable 
is  based  upon  the  fact  that  if  we  toss  the  coin  many  times  then  roughly 
half  of  the  tosses  yield  heads  and  half  yield  tails.  In  addition,  there 
are  two  analogous  theorecical  results  proved  or.  the  basis  of  our  model, 
the  weak  law  of  large  numbers  and  the  strong  la>.  of  large  numbers  which 
strengthens  the  bond  between  intuition,  experiment,  and  theory. 


The  cerumen  model  for  two  independent  tosses  of  a  fair  coin  is 
almost  as  simple.  The  probability  space  can  be  taken  to  be  as  follows: 

(l  =  -{HH,  HT.  TH,  TT}  ,  ¥  consists. of  all  subsets  of  and  if 
Aefi("e"  meaning  "an  element  of")  then  p(A)  =  (Number  of  elements  of  A) A. 
Suppose  we  ask  now  "what  is  the  probability  of  obtaining  both  a  head  and 
a  tail?"  In  our  model  this  means  "what  subset  A  corresponds  to  'both  a 
head  and  a  tail'  and  what  is  P(A)?"  The  answer  is  A  =  -[HT,  TH}  and  P(A) 

=  2A  =  1/2. 

Similar  examples  can  he  given  for  dealing  cards,  drawing 
marbles  out  of  urns  and  many  other  real  life  situations  where  there  is 
an  element  of  chance  involved.  We  refer  the  reader  to  a  book  on 
probability  theory  such  as  (2,  Vol.l)  if  his  interest  moves  him  to  find 
out  more  about  them.  It  is  important  for  an  understanding  of  the  subject 
to  discuss  one  example  which  is  essentially  different  from  the  coin 
tossing  experiment. 


2.1.3  Dart  Throwing 

The  coin  tossing  example  had  the  property  that  P(A)  >  0 
whenever  A  contained  only  a  single  element  of  ft.  If  this  were  always 
true  it  would  be  unnecessary  to  introduce  the  collection  ¥  into  the 
model.  We  now  describe  a  situation  where  it  isn't.  Suppose  a  dart  is 
tossed  at  an  infinite  dartboard  equipped  with  a.  cartesian  coordinate 
system.  It  is  aimed  at  the  origin.  As  any  enterprising  dart  thrower 
realizes,  the  chances  of  hitting  exactly  the  origin,  or  exactly  any 
otb^r  point  is  zero.  On  the  other  hand,  the  probability  of  hitting  a 
nonzero  area  on  the  dartboard  is  general.! y  non  zero. 

For  this  example  let  ft  cons-5  at  of  all  the  points  on  the 
dartboard,  or  better  yet,  all  pairs  of  real  numbers  (x,y)  which  will 
he  coordinates  of  points  on  the  dartboard  (bull's  eye  at  the  origin). 

Let  t  be  the  smallest  a-field  of  sets  containing  all  rectangles.  The 
satisfying  this  description  turns  out  to  be  quite  large  and  contains, 
for  example,  all  regions  whose  boundaries  are  continuous  curves.  Finally 
we  let 

.P(A)  ]V(x2+y2)/c2dxdy 

A 

Then  P (<p)  =  0,  P(ft)  =  1,  and  P  satisfies  the  additivity  condition  because 
of  the  properties  of  the  integral.  Furthermore,  p(A)  =  0  whenever  A  is 
a  single  point  in  the  plane  or,  for  that  matter,  even  a  smooth  curve. 

The  function  is  called  the  probability  density  or  distribution  density 
function.  The  constant  c  is  related  to  the  skill  of  the  dart  thrower 
(the  smaller  c  the  better).  It  has  been  found  that  this  particular  type 
of  density  function  describes  an  actual  dart  throwing  experiment  part¬ 
icularly  well.  That  is,  if  the  dart  is  thrown  N  times,  where  N  is  large, 
and  the  number  of  times  a  hit  is  made  inside  a  given  set  A  is  NA,  then 
N^/N  approximates  p(A)  and  the  approximation  improves  as  H  increases. 


2.1.4  Conditional  Probability  and  the  Urn .Model 


The  machinery  of  sample  spaces  is  ideally  suited  to  deal 
with  problems  involving  what  is  known  as  conditional  probability.  Lot 
us  consider  a  problem  where  this  concept  is  involved.  Suppose  we  have 
two  urns  numbered  1  and  2  and  each  contains  a  red  ball  and  a  black  ball. 

A  ball  is  drawn  from  urn  1  and  placed  (without  looking  at  the  colour) 
into  urn  2.  Then  a  ball  is  drawn  from  urn  2  and. it  is  black.  What  is 
the  probability  that  the  ball  originally  drawn  from  urn  1  was  black? 

To  solve  this  problem  notice  first  that  there  are  two  possible  compps- 
itions  for  urn  2,  each  equally  likely  R^I^and  B1R2B2  (to  clarify 
matters  we  have  used  subscripts  to  denote  the  urn  of  origin  of  the  various 
balls).  Let  ft  =  ^Ri ,R2jb2»b1 >^2>B2  1  j  ¥  consist  of  all  subsets  of  ft  and 
P  be  defined  by  p(A)-  =  (number  of  elements  in  A)/6.  This  assignment  of 
probabilities  is  consistent  with  our  intuitive  prejudice  that  each  ball 
in  a  given  urn  has  an  equal  chance  of  being  drawn.  The  subset  of  ft 
corresponding  to  the  event  "a  black  ball  is  drawn  from  the  second  urn" 
is  A  =•  {B2,Bf,B2}  .  Since  we  know  this  event  has  occurred  we  restrict 
our  attention  to  a  new  sample  space  (ft,  i,Pi)  where  ftj  =  A,  ¥1  consists 
of  all  subsets  of  ft*  and,  since  each  single  element  set  has  equal 
probability  in  (ft,^,P),  we  ask  for  the  same  property  in  (fti ,  ^  1  ,P'i ) . 

This  means  Pi(B)  =  (ft  elements  in  B)/3  whenever  Be  \ .  Such  a 
definition  leads  to  the  conclusion  that  the  probability  of  a  black  ball 
on  the  first  draw,  which  corresponds  to  the  primed  letters  in  A,  is 
2/3.  Experience  shows  2/3  is  correct,  in  support  of  the  usefulness  of 
this  definition  of  sample  space. 

In  general  if  (ft,^,P)  is  a  sample  space  and  ft*,  with 
P(ftl)  >  0  is  a  subset  specifying  a  condition  known  to  be  satisfied, 
then  the  conditional  sample  space  is  defined  to  be  (fti , 1 ,  Pi)  where 
=  (ftnA  :  Ae*  }  and  Pj(B)  =  P(B)/P(A)  when  B  e  ¥j,. 


2.2  Random  Variables 

2.2.1  Definition  and  Examples 

Having  described  what  is  meant  by  a  sample  space,  we  are  now 
ready  to  discuss  random  variables,  a  concept  vrJLch  is  indispensible  in 
the  study  of  Monte  Carlo  simulators.  A  random  variable  X  on  a  prob¬ 
ability  space  (ft,  "4  ,P)  is  a  function  whose  domain  is  the  set  ft  and 
whose  range  is  in  the  real  numbers  such  that  (10  :  a  <  X(o>)  <  8}  e  ^ 
for  all  real  numbers  a  <  8.  In  practice  it  is  useful  to  deal  with  random 
variables  whose  range  is  an  n-dimensional  Euclidian  space,  but  we  shall 
avoid  this  complication  for  the  moment.  The  probability  that  a  random 
variable  assumes  a  value  in  an  interval  (a, 3)  is  denoted  by  P(a  <  X  <  6) 
and  given  by  P(A)  where  A  is  the  set  :  a  <  X  (w)  <  8^  ,  i.e.  the  set 
of  elements  of  ft  which  lead  to  a  value  of  X  between  a  and  8.  It  is  clear 
now  why  we  assumed  such  sets  were  in  ¥  since  is  the  domain  of  the 
probability  function  P, 

"Natural1'  examples  of  random  variables  abound  in  real  exper¬ 
iments  where  chance  is  involved.  If  a  coin  is  tossed  N  times  the  total 
number  of  heads  is  a  r.indom  variable.  In  the  dart  throwing  experiment 
the  x  coordinate  of  a  hit  is  a  random  variable,  so  is  the  y  coordinate 
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T/P 

and  the  distance  from  the  origin  (x2  +  y2)  7  .  There  will  be  many  more 
examples  in  the  Monte  Carlo  methods  we  are  about  to  discuss. 


2,2,2  Distribution  and  Density  Functions 

In  working  with  a  single .random  variable  X  it  is  occasionally 
useful  to  devise  a  new  sample  space  (ft  ,  ^  , -P  )  which  is  simpler  in 
many  cases  than  the  original  one. and  yft  leads  £o  the  same  conclusions 
as  far  as  probability  theory  is  concerned,  This  sample  space  is  defined 
as  follows :  ft  is  the  real  numbers  5  ¥  consists  of  all  the  Borel- 
measurable  subsets  of  the  real  numbers  and  P  is  completely  specified 
by  giving  its  value  on  every  interval  (a,  8): 

Px((a,  e))'=  P(a  <  X  <  8) 

Px  called  the  probability  distribution  function  for  X. 

In  many  cases  it  is  possible  to  find  a  real  valued  function 
1 fix )  with  the  property  that 

Px((ot,  8))  =  j  <?(x)  dx 
a 

for  all  (a,  8).  Then  ^  is  called  the  probability  density  function  for 
the  random  variable  X. 


2.3  Expectation  and  Variance 

We  now  come  to  an  exceedingly  useful  concept.  The  expectation 
of  a  random  variable  corresponds  to  the  answer  wo  obtain  if  we  make 
many  observations  of  the  random  variable  and  take  their  average.  It 
is  defined  most  generally  as  follows: 

E(x)  =  j  X(w)  P(dw) , 
ft 

The  expression  on  the  right  is  an  integral  over  the  probability  space 
and  is  defined  to  be  the  limit  of  finite  sums  of  the  form 

?  \  P(V 

k 

where  ,  is  a  disjoint  decomposition  of  ft  and  the  s^  are  chosen  in  such 
a  way  that  the  step  function  with  constant  value  s^  on  (k  =  1,  2,,.,) 
approximates  X  in  a  certain  sense.  The  theory  we  ‘nave  <lealt  with  so 
far  is  not  equal  to  the  task  of  making  this  definition  precise. 

We  shall  therefore  give  another  definition  which  is  equivalent 
to  the  above  one  but  applies  only  in  the  case  where  the  random  variable 
X  possesses  a  density  function  (p  .  Then 


15 


E(X)  =  i  x<P(x)  ax. 

'  _ m 


One  of  the  annoying  faets  of  the  business  is  that  the  integral  in  either 
definition  may  fail  to  exist.  There  is  nothing  that  can  he  done  about 
this ,  and  ve  shali  not  deal  with  random  variables  whose .  expectation 
fails  to  exist  in  the  sequel.  Expectations  and  random  variables  have  a 
number -of  useful  algebraic  properties;  we  now  state  a  most  basic  one. 

Let  Xj  ‘and  X2  be  any  two  random  variables  whose  expectation  exists  and 
cf  be  a  real  number.  Then  aXj  +  X'2  is  again  a  random  variable  whose 
expectation  exists  and 

E(aXi  +  X2)  =  ciE(Xj)  +  E(X2) 

The  expectation  E(X)  is  often  called  the  mean  of  X.  This  is 
because  if  we  perform  an  experiment  whose  outcome  is  X  N  times  in  such 
a  way  that  each  repetition  is  independent  of  the  previous  ones  and  take 
the  average  of  the  values  of  X  we  obtain  in  this  way,  this  average  turns 
out  to  get  closer  and  closer  to  E(X)  as  N  gets  large.  This  process  of 
repeating  an  experiment  is  amenable  to  a  much  more  careful  analysis 
which  we  postpone  momentarily. 

Instead  we  pursue  another  idea.  Just  as  it  is  useful  to  know 
the  mean  of  a  random  variable  X,  it  is  useful  to  have  a  measure  of  how 
much  the  values  of  X  can  he  expected  to  deviate  from  the  mean.  The 
concepts  which  have  been  introduced  for  this  purpose  are  the  variance 
and  the  standard  deviation  of  X.  The  variance  is  defined  to  be  the 
expectation  of  the  random  variable 


(X  -  E(x))2  a  X2  -  2  X  E(X)  +  (E(X))2 


That  is, 


Var  (X)  =  E(X2)  -  2  (E{X))2  +  (E(X))2  =  E(X2)  t  (-E(X))2 
In  case  X  has  a  probability  density  function  jb  (x) 

Var  (X)  -  |  (x  -  E(x))2  (x)  dx 


The  standard  deviation  of  X  is  defined  to  be  the  square  root  of  the 
variance,  i.e. 

0  (X)  =  (Var  (X))1/2 

In  order  to  make  our  statement  that  the  variance  and  standard  deviation 
gv-  -5  us  an  idea  of  how  much  the  values  of  X  can  be  expected  to  deviate 
from  the  mean,  we  state  a  result  known  as  the  Chebychev  inequality 
[l,  Prop.  1.7,  p.  4]  : 

P  ([  10  :  |y'w)  ~  E  ( X ) | -/ a  >  e}  )  <  e~2 
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2.4  Independence  r 

2.4.1  Product . Spaces 

In-order  to  "build  a  model  to  handle  such  situations  as 
repeating  a  chance  experiment  we  introduce  the  concept  of  product  of 
independent  probability  spaces.  To  this  end  let  -P^)  k  =  1, 

...,n  be  probability  spaces.  Let  Q  be  the  sample  space  whose  elements 
are  n-tuples  of  elements  of  k  =  N.  That  is, 

=  1,...,N} 

or  equi valently 

fi  =  %  x  ft2  x  ...  X 

We  let  ¥  be  the  smallest  a-field  containing  all  sets  of  the  form 
A  =  A1x  A2x  ...x  Ajj, 


where 

\  e  ^k’  k  = 

We  define  P  first  on  sets  of  the  form  A  =  Aj  x  . . .  x  A^  by  setting 


P(A)  =  Pi  (AO  P2(A2)  ....  P^) 


and  then  extend  it  to  all  of  ¥  by  making  us?  of  the  addition  axiom 
and  of  continuity  properties  of  the  P.  .  It  can  be  checked  that 
(W,  ,  P)  is  a  probability  space.  We  shall  call  it  the  .product  of 

the  probability  spaces  (&.,  ¥•  ,  P.  ),  k  =  1,...,  N,  and  denote  it 
occasionally  by 


N 

n 

k=l 


<V 


V- 


We  shall  also  have  occasion  to  make  use  of  an  infinite  product 


n 

k=l 


of  sample  spaces.  The  concept  and  definition  are  analogous  but .too 
complicated  to  attempt  here. 


2,4,2  Motivation  and  Definition 

t 

of  the  form  fij  x  ...  x  fi  x  A^xfi  x  ...  x  which 
exhibit  a  very  important "property  -  namely, 


Sets 

A 

we  denote  by  A^ 


ny  “p(v  -p(y 
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whenever  kj  #  -k2  #  • .  •  -  t  k£ .  -It  is  important  enough  to  have  a  definition 
all  its  own.  Events  Bj,  32,.  ...  B.  are  said  to  he  independent  of 

* 

?  {Bjfl  ...  HB£)  =  P  (Bi)  P  (B2)  ...  P  {3,) 

The  notion  of  independence  in  our  model-  is  a  rigorous  statement  of  the 
notion  as  derived  from  experiments  involving  chance.  We  shall  make  much 
use  of  its  consequences  for  random  variables,  end  we  proceed  to  derive 
some  of  them  presently. 


If  Xj  is  a  random  variable  on  (fij,  'f  j,  Pi),  then  it  can  he 
treated  as  a  random  variable  on 


(«,?,  p,)  =  n  (a,*.,.?) 

k=l  *  K  5 


a  natural  way  -by  letting 


...»  -  Xi(o)i) 

We  find  then  that 

P  (a  <  X!  <  8)  =  ?  («  <  Xx  <  g) 

so  that  as  far  as  probabilities  are  concerned  Xj  and^  are  identical. 
Suppose  3^  is  a  random  variable  on  (fl.  , ¥  ,  P  )  and  X,  is  the  corresponding 
random  variable  on  the  product  space,  k  -1,  ...  ,  N.  Then 


PCX^  eJ1>&...&Xk  eJt)sP(^eJ'1)...P  s  J£) 

where  kj  ?  k2  f  . . .  f  k£  and  J. ,  i  =  1,  . . . ,  SL  are  Borel  sets  of  real 
numbers.  This  follows  easily  x'rojn  the  independence  property  of  events 


which  we  discussed  in  the  preceding  paragraph.  Hence  we  define  the 


notion  of  independence  for  any  £  random  variables  Xj ,  . . . ,  X  defined 
on  a  common  probability  space  by  the  statement:  " 


l 

P  (Xi  e  J]  &,  ...,  &  X  e  J  )  =  H  P  (X.  e  J.). 

*  *  i=l  11 


2.4.3  Independence  and  Expectation 


they  are  in 


If  Xi ,  . . . ,  Xj.  are  random  variables ,  then  so  is  il  X.  .  If 
in  addition  independent  it  can  be  proved  that  k=l 


1  independent  it  can  1 
N  N 

E  (n  X.  )'=  n  E  (X,  ). 

1 _ T  1 _ -I  ^  • 


For  our  purposes  this  property  is  not  as  important  as  one  of  its  consequences 
—  namely 

N  N 

Var  (£  \ )  =  l  Var 
k=l  *  k=l 
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■3^2.02  Z 


which  we  prove  here  wham  2  =  2  for  tbs  sake  of  iH 

Vsr  (Sj  *  S-)  =  BS(Si  +  %z )2)  -  ls{?i  +  r2>32 
=  SCS-.2  -5-  2Sj  22  +  %2)  "  (CS(Xi))2  *  2E(%)  2(>:2)  -r  (S{a2))2} 

-  E{X|2)  +  £S(Xj  t2)  4-  S{X22}  -  Cs(xs)2  -  22(Si)  2(S2)  -  Cs{52))2 
=  s(s22)  -  Cs(x2))2  -  s(x22)  -  (b(s2))2  =  ¥22-  (:q)  *  ?a?  (x2). 

2.5.  The  Uniform  and  Sormal  Bjstribiitlcns 

Sc  far-  *«e  cave  dealt  with  a  large  susbsr  of  definitions  sue. 
theory.  In  order  to  rake  use  of  this  theory  we  seed  to  have  sees  sore 
specific  knowledge  about  probability  spaces  which  recur  frequently  in 
building  statistical  models  for  real  situations,  We  shall  describe 
wo  such  spaces  now. 


Tbs  first  of  these  is  known  as  the  uniform  distribution  and 
is  defined  as  follows:  let  S  be  the  interval  (a,s>)  .and  let  ¥  consist 
of  all  the  Borel  measurable  subsets  of  (a.b).  The  probability  function 
is  defined  to  be 


?(A)  = 


l 


dx 

b-a 


where  — is  the-  density  function.  We  then  calculate  the  mean 


and  the  variance 


and  remember  these  facts  for  future  reference. 


The  second  space  we  must  mention  is  the  normal  distribution 
which  is  defined  as  follows:  Let  ft  be  the  whole  real  line  and 

let  ¥  consist  of  all  the  Borel  measurable  subsets  of  the  real  line. 
We  fix  values  of  the  parameters  v  and  a  and  define  the  probability  of 
a  set  A  by 

so  that  the  density  function  is 

1  ~{x-y)2/2a2 

e 
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a 


The  mean  them  is  gives  by 

212?  =  J" 

*g|j®  V52T2SnCS  "fcy 


-(z-is)2/232 


iffy 


*  _ SB 


fix  =  o2 


froa  which  the  stsr.fi.erd  deviation  is  a-  The  preceding  integrals  cannot 
fee  evaluated  isj  eiesesSary  seans  and  one  usually  trassfcms  then  into 
linear  combinations  of 


which  are  ccnmonly  found  is  tables  of  integrals. 


2.6-  The  Central  Licit  Theorem 

This  pretty  well  completes  the  outline  of  the  sost  essential 
statistical  ideas  necessary  to  plunge  into  a  study  of  r-fente  Carlo  methods. 
There  remains  unsent icned,  however,  a  fundamental  and  far  reaching 
theorem  which  proves  invaluable  in  understanding  the  answers  a  tonte 
Carlo  calculation  cranks  out.  The  central  limit  theorem  assumes  we 
have  an  unlimited  number  of  independent  random  variables  X  (k=l,2, . . . ) 
and  each  has  the  same  distribution;  that  is,  given  (a,S),  then 

p(a  <  X.  <  6)  is  the  same  for  all  k=l,2, _ leading  to  a  mean  y  and 

a  standard  deviation  a-  Then  we  can  conclude  that  for  If  sufficiently 
large,  the  random  variable 


1 

N 


H 

Z 

k=l 


*k 


has  an  approximately  normal  distribution  with  mean  y  and  standard 
deviation 


o 

Jn 

Strictly  speaking,  the  terj  "approximately"  above  requires  more  precise 
definition.  We  therefore  refer  the  reader  to  [1,  p.  170]  where  the 
following  precise  statement  nan  he  found: 

Theorem 


Let  Xl5  Xj ,  ...  be  independent,  identically  distributed  random  variables 
each  of  mean  y  and  variance  o.  Then 
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s 

n 


Ij  *  %2  *  ---  X 


a  JE 


S> 

■*  1  (0,1) 


Tbs  symbol  ®  S  (0,1)  stasis  for  converges  in  distribution 
to  tbs  norm si  distribution  of  ^<g~r»  zero  standard,  deviation  1.  It 
seams  that  if  Fn  (a)  =  FfSj,  <  a)  is  the  distribution  function  of  %  then 


lin 

j>»o 


2-7  Sistograns 


Suppose  vs  sre  given  a  device  vhich  is  able  to  generate 
randoms  numbers  chosen  fron  a  distribution  vith  Q  =  (a,b)  and  density 
^pfe)-  In  this  section  ve  discuss  tee  problem  of  obtaining  an  approx¬ 
imation  to  c?(x)  by  platting  a  histogram  iron  number  S?  of  random 
niitbers  (r )  produced  by  our  generator. 


First  ve  state  vhat  is  react  by  a  histogram.  Let  a  partition 
{a  =  x  <  xi  <  . « .  <  x  =  d}  of  the  interval  (s,b)  be  fixed-  Thus 
is  subdivided  into  cells  a^  =  (x^_i,  xjj)  k  =  1,...,  n.  The  approx¬ 
imation  ^«tx)  is  defined  to  be  the  step  function  vhose  constant  value 
on  the  interval  £lv  is  obtained  by  counting  the  number  of  the  r^  vhich 
fall  into  this  interval  and  dividing  the  result  by  If(x^  -  x^  1 ) - 


There  is  another  way  of  describing  tba  histogram  vhich  is 
core  useful  for  theoretical  purposes.  Let  x,  be  the  characteristic 
function  of  the  cell  a,  ,  that  is, 

k  V 


Xk(x) 


1 

0 


X  E  SI, 
k 

otherwise 


V. 


Then  is  a  random  variable  on  (a,  V* ,  P).  We  can  imagine  to  be 
a  random  variable  on  the  i-th  factor  of  the  product  space 

N 

n  (n,*,  p) 

l 

and  denote  its  extension  to  the  whole  product-  space  by  x^  i*  For 
fixed  k  then  -  i  =  1,...^  N  are  identically  distributed  independent 
random  variable^.  Furthermore,  if  we  let  <jf>  denote  the  value  of 
^(x)  on  0^,  we  obtain  ’ 


1 

Vs’k  =  «VV?  A 

^  k  is  also  a  random  variable  on  the  product  space.  We  calculate 
the ’mean  and  variance  of  y.  : 

it 


T 
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H(X_)  =  f  1 SPJx)  ex  =  V^i 
¥sr  (\)  =  J  I2  P  (x)  cz  -  y^2  =  yjl  -  Pj.) 


It  follows  them  that 


=r  v>  =  Wl)2 


and 


°'£  H,k:  ■ 


‘“k 


? — ) 
X-l' 


\ 

If  H  is  large  then  ve  can  use*  the  central  Unit  theorem  to  conclude 
that/Sg -u:  vill  he  approximately  normally  distrihuted  implying  that 
it  vill  fall  within  two  standard  deviations  of  about 

98/>  of  the  tine.  It  is  perhaps  more  illuminating  to  look  at  the 
"relative  deviation" 


“  WV^2 


Then  if  pv«l ,  which  must  he  the  case  if  the  nartition 
points  are  finely  spaced  in  order  to  make  the  histogram  show  detail, 
the  relative  deviation  is  approximately  or  (expected  number 

of  counts  in  cell  This  is  a  useful  rule  of  thumb  to  remember 

when  designing  Monte  Carlo  algorithms.  It  is  heavily  dependent,  of 
course,  on  the  assumptions  of  independence  and  small  %. 

We  have  worked  out  the  histogram  theory  for  a  distribution 
with  fi  =  (a,b)  and  probability  density  .  The  theory  is  not  really 
restricted  to  these  assumptions.  More  generally  we  can  start  from 
an  arbitrary  probability  space  p)  and  choose  a  partition 


n  = 


N 


k=l 


ftk  disjoint  and  fij-E-f,  k=l,2,...n  of  the  set  fi.  Again  we  choose  a 
sample  (uj,  102 ,  •  •  •  sulj)  of  N  random  uoints  from  the  distribution 


Since  P  may  not  be  given  by  a  density  function  it  is  not 
possible  to  define  the  histogram  as  a  stepfunction  approximation  to 
the  density  function.  Instead  we  let  pn  ^  equal  the  fraction  of 
the  falling  into  U k  and  consider  it  an  approximation  to 
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-lai?  =  V 

We  define  as  cefere  ana  find  that 

E 

=  S  Xk,z 

is  a  random  variable  on  tee  space 

B 

5  P) 

ifiiose  T^n  turns  osu  to  "be  ana.  stasdard  deviation  turns  out  no  "be 

Vk(l-pk)H"1/2 

Again  the  relative  standard  deviation  is 

(l/pk-l)1/2lf1/2 

and  our  rule  of  nhurfo  that  relative  standard  deviation  ~  (expected 

number  of  counts  in  ft.  )“-/2  still  auulies. 

s. 

2.8  Algorithms  and  Monte  Carlo 
2.8.1  General  Structure 


We  shall  not  define  precisely  what  is  meant  "by  an  algorithm. 
It  is  sufficient  for  our  purposes  to  say  that  an  algorithm  is  a  finite 
sequence  of  operations  and  of  things  to  he  operated  on,  instructions 
giving,  the  order  in  which  the  operations  are  to  be  performed  and 
when  the  execution  is  to  stop,  together  with  a  guarantee  that  the 
operations  can  always  be  carried  out  under  stated  conditions.  A 
computer  program,  of  course,  fits  the  description  we  have  just  given 
and  it  is  in  this  context  that  we  shall  be  discussing  algorithms  from 
now  on. 


Since  we  shall  not  be  discussing  any  particular  algorithm 
in  the  next  few  paragraphs  and  consequently  shall  not  be  concerned 
with  the  internal  details  of  how  a  particular  algorithm  is  put  together, 
it  is. appropriate  to  use  the  well  known  function  notation  as  a 
shorthand  for  an  arbitrary  given  algorithm.  That  is,  let  U  be  the 
vector  of  numbers  which  the  algorithm  uses  as  initial  parameters  and 
V  be  a  vector  whose  components  are  the  numbers  which  it  produces.  Then 
we  write 


V  =  $  (U) 

to  denote  the  action  of  the  algorithm. 
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Generally  algorithms  are  demised  so  tbst  ve  can  vise  a 
oimaaatar  ao  calculate  some  function  which  is  at  least  theoretically 
hmiB.  ?sr  example,  let  A  be  an  n  x  n  matrix.  of  complex  numbers. . 
‘Then  it  is  known  that  A  has  a  complex  eigenvalues.  Alternatively 
ve  could  say  that  there  exists  a  function  F\A)  whose  value  is  a 
vector  cf  a  components  —  the  eigenvalues  of  A  in  sene  order.  It 
aa  an.  important  problem  to  devise  an  algorithm  $(A,k)  which  produces 
a  close  approximation  to  the  eigenvalues  of  A.  Ve  ash  that  the 
algorithm  at  least  satisfy 


^  $(A,k)  =  F(A) 

Experience  tells  us  that  such  an  algorithm  would  utilize  seme 
iterative  method  and  that  the  time  T(k)  required  to  execute  the 
algorithm  would  go  to  =>  as  k  -*»,  ‘That  is  something  one  must  live 
with. 


With  this  preamble  ve  now  turn  our  attention  to  what  we 
shall  call  a  Monte  Carlo  algorithm.  We  describe  it  at  first  in  a 
somewhat  idealized  form  in  order  to  clarify  what  is  meant. 


Suppose  that  our  computer  has  attached  to  it  a  source  of 
random  numbers  whose  distribution  is  known  so  that  in  addition  to 
the  usual  slate  of  instructions  there  is  one  which  enables  us  to 
introduce  a  number  from  the  random  source  into  the  computation.  If 
an  algorithm  makes  use  of  these  random  numbers,  we  should  intuitively 
expect  that  the  answer  will  also  be  described  by  a  probability  distr¬ 
ibution.  Ve  can  describe  a  Monte  Carlo  algorithm  in  terms  of  our 
function  notation  if  we  imagine  that  instead  of  having  a  random 
number  generator  built  into  the  computer  we  allow  the  computer  to 
have  access  to  any  finite  part  of  an  infinite  sequence  of  numbers 
produced  separately  by  a  random  number  generator.  Let  such  a 
sequence  be  denoted  by  r  =  (rj,  r2-.-).  The  words  "finite  part" 
in  the  second  to  last  sentence  deserve  emphasis.  We  wish  to  exclude 
such  functions  as 


"Compute  the  mean  of  the  sequence  r" 


or 


"If  lim  r,  is  greater  than  10  commence  step  3". 

Fortunately  in  most  algorithms  it  is  possible  to  place  an  upper  bound 
N  on  the  length  of  the  sequence  required  by  the  algorithm.  Then  it 
suffices  to  work  with  r  =  (rj,  r2,...r^). 

For  theoretical  purposes  we  let  (ft  ,  f  ,  P  )  be  the  prob¬ 
ability  space  describing  the  random  number  generator?  Then  r  may  be 
treated  as  a  draw  from  the  probability  space 

N 

(«,¥,P)  =  Jl  (fl0>VP0)’ 

1 


or 


n 

4=1 


2k 


as  the  case  nay  be-  The  Monte  Carlo  algorithm  Is  denoted  by 


V  =  $(U,?,k). 


The  parameter  k  ba.a  been  introduced  to  allov  for  convergence 
properties.  It  -is  evident  that  for  fixed  U,  V  is  a  vector  valued 
random  variable  on  (fi,  ,P) .  It  is  therefore  reasonable  tc  ask  for 
the  probability  distribution  of  V.  Generally,  tnis  is  too  difficult 
a  question  to  answer,  and  for  practical  purposes  certain  properties 
of  this  distribution  suffice. 


2.8;2  Monte  Carlo  Integration 

It  is  time  now  to  describe  an  algorithm  which  is,  historically 
speaking, . a. prototype  for  all  Monte  Carlo  calculations.  The  purpose 
of  this .algorithm  will  be  to  approximate  the  integral  of  a  given 
function  f(x)  on  the  interval  (0,1).  It  is  presumed  that  f  is  given 
to  us  try  specifying  d  parameters  -  they  could  be  the  coefficients  of 
a  polynomial  of  order  d  -  1,  or  they  might  be  the  values  at  d  points 
in  the  interval  (C,l)  if  f  is  assumed  precisewise  linear.  The  exact 
method  of  specification  is  of  no  direct  concern  to  us.  Let  us  simply 
assume  that  once  the  d  parameters  are  given,  we  know  how  to  calculate 
f(x)  for  any  x  e-(0,l). 


The  actual  computation  of  the  integral  proceeds  as  follows. 
Given  N,  choose  fl  independent  random  numbers  r  =  (rj,  r2, . . . jr^)  from 
the  uniform  distribution  on  (0,1)  and  then  compute 


1  N 

I  =  $  (f»r,N)  =  ±  E  f(r  ) 
k=l 


we  claim  that  I  in  an  approximation  to 


dx 


in  the  following  sense: 


E(I)  =  f  E  E(f) 
k=l 


f(x)  dx  =  vi 


Var  (I)  =  | 
a  (i)  =  if 1/2 


We  summarize:  The  mean  of  $(fjr,N)  is  the  required  integral  and  its 
standard  deviation  tends  to  zero  as  Furthermore,  $  is  an  average 

of  random  variables  that  for  sufficiently  large  N  we  can  assume  it 
will  be  approximately  normally  distributed  with  mean  y  and  standard 

deviation  C  in  particular  the  probability  that 

\ 

j  $  -  y  |  <  2  C  if 1/2 
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is  -about.  0.98.  frcm^  which  we  can  conclude  that  with  probability  0.98 
will  approximate  f~  fix')  dx  to  within  e  it  S  >  (2 C/s)2. 

0 

Before  leaving  this  example,  it  is  import -art  to  notice 
that .tbe .number  of  computations  of  f  and  hence  the  time  taken  to 
execute. the  algorithm  increases  as  I?  hut  the  standard  deviation 
only. decreases  as  2i-l/S.  If  we  compare  this  with  a  quadratically 
convergent  standard  algorithm,  ve-can  see  that  Monte  Carlo  is  a 
preposterously .expensive  way  to  get  accuracy. 

We  can  also  use  this  example  as  a  guide  to  formulate  a 
criterion  for  convergence  of  a  Monte  Carlo  algorithm.  We  say  that 
the  algorithm 


$  (U,  r,  n) 


converges  to  F(U)  as  n-*=  if 

Pn  =  E(*(U,  *  ,n) )  -*■  F(U) 


and 


Var  ($(U,-,n))  0. 

Unfortunately,  neither  this  criterion,  nor  any  alternative 
has  been  proved  for  the  Monte  Carlo  algorithms  proposed  to  date  in 
Kinetic  Theory.  We  state  it  here  primarily  to  throw  light  on  the 
nature  of  the  answer  a  Monte  Carlo  algorithm  produces. 


2.9  Random  Number  Generators 


2.9.1  The  Uniform  Distribution 

We  stated  at  the  beginning  of  the  preceding  section  on 
Monte  Carlo  algorithms  that  the  random  number  generator  attached 
to  our  computer  was  an  idealization.  In  practice  a  substitute  is 
used. 


There  are  several  standard  algorithms  known  which  have  the 
property  that  they  produce  a  sequence  of  numbers  which  as  far  as  its 
statistical  properties  are  concerned  could  have  been  produced  by 
sampling. from. the  uniform  distribution  on  (o,l)  (there  is  nothing 
basic  about  this  interval,  it  has  been  adopted  by  convention  mainly). 
The  algorithms  must  be  given  a  starting  value  r  so  that  we  could 
write  a  typical  one  as 


R(r  n) 

0, 

which  gives  the  n-th  pseudo  random  number  if  the  starting  number  is 
r  .  Among  the  statistical  properties,  we  demand  that  histograms 
built  from  the  sequence  R(r  ,n)  be  statistically  indistinguishable 
from  a  histogram  built  from0 the  same  number  of  samples  from  a  uniform 
distribution. 


i 
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It  is  tempting,  "but  inconsistent  with  our  goals,  to  go  into 
detailed  descriptions  of  the  tests  which  are  used  in  this  connection. 

We  refer  the  reader  instead  to  [3,  Ch.12] . 

The  test  methods  known  to  date  for  generating  uniformly 
distributed  random  numbers  are  known  as  the  recursive  congruential 
generators  (5).  The  n-th  pseudo-random  number  R(rQ,n)  is  given  by 

br  n  =  R(r . ,n)  (modulo  m) 
o  o 

for  a  suitable  choice  of  the. parameters  (b,m).  This  choice  requires  a 
considerable  amount  of  number  theory  to  make  properly.  The  congruential 
methods  are  both  fast  and  pass  the  statistical  tests  with  flying  colours. 

There .are  two  basic  methods  for  devising  an  algorithm  to 
generate  random  numbers  from  a  distribution  given  by  an  arbitrary 
density  function  jb(x)  on  which  make  use  of  a  source  of  uniformly 

distributed  random  numbers. 


2.9.2  Method  I 
Let 

*(*)  «  f 

J  -CO 

$  is  a  nondecreasing  continuous  function  which  assumes  every  value 
in  the  interval  (0,1).  Consequently,  there  exists  at  least  one 
nondecreasing  inverse  function  on  the  domain  (0,1)  such  that 

$($-1(y))  =  y. 

<i>"1  is  a  random  variable  on  the  uniform  probability  space  with 
fi  =  (0,1).  Since 


xi  <_  <J>-1(y)  £  x2 


implies 


^(xl)  £y  i  *(x2) 


it  follows  that 

P(xx  <_  <J>_1  <_  x2)  =  P( 4>(xi )  <  y  <  $(x2))  -  4>(x2)  -  4>(xi ) 

Hence  $  1  is  distributed  according  to  the  density  function ^(x) ,  and 
if  R(r0,n)  is  a  random  number  generator  for  the  uniform  distribution 
on  (0,1)  then 

i~1('R(rQin)) 

is  a  random  number  generator  for  the  distribution  with  density £ (x) . 
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Method  I,  though  completely  general  in  theory,  is  :not  always 
useful  from  a  practical  point  of  view.  The-  problem  arises  when  there 
is  no. easily. computable  expression  for  $  1  or,  even  worse,  if  we  are 
given  <p  and -there  is  no  easily  computable  expression  for  either  $  or 
$  1.  In  such  cases  Method  II  may  be  preferable. 


2.9.3  Method  II 

Let  p  be  a  given  function  not  identically  zero  on  the  interval 
(a,b)  and  satisfying  the  inequality 

0  <_  f{x)  <_  c. 

We  shall  describe  a  meiwu.  for  obtaining  a  random  number  from  the 
distribution  defined  by  the  'density 

<5p(x)/  I  f>(x)  dx 
^  a 

Let  x  be  a  random  number  from  the  uniform  distribution  on 
(a,b)  and  y  be "a  random  number  drawn  from  the  uniform  distribution  on 
(0,c).  Then 

r* 

P(  <f(x)  1  y)  =  I  ?(x)  dx/(b-a)c 
J  a 

and  the  distribution  for  x  given  that  the  condition  <p(x)  <  y  is 
satisfied  is 


P(a  <  x  <  8  on  the  condition  ^(x)  <  y) 


That  is,  the  above  conditional  distribution  for  x  is  the  distribution 
we  seek  to  compute. 


On  the-  basis  of  the  theory  just  derived,  we  devise  the  random 
number  algorithm  as  follows: 

A.  Choose  rj  and  r£  from  t’  miform  distribution  on  (0,1).  Then 

x  =  a  +  (b-a)  rj ,  and  y  =  cr£ 

will  be  uniformly  distributed  on  (a,b)  and  (0,c)  respectively. 

B.  If  ^(x)  <_  ythen  x  is  the  required  random  number  and  ye  are  finished, 
Otherwise  return  to  A  to  repeat  the  process. 


One  of  the  interesting  features  of  this  algorithm  is  that  the 
number  of  steps  necessary  to  generate  a  random  number  x  is  itself  a 
random  variable.  Some  further  analysis  is  in  order.  The  probability 
that  x  is  accepted  on  the  first  pass  is 


.xjc 


dx/c(b-a). 
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The  probability. of . acceptance  on  exactly  the  n-th  pass  is  the  same 
as  the  probability  of  -non-acceptance  on  the  first  (n-1)  passes  and 
acceptance  on  the  n-th  pass,  or 

(l-p)n-1p. 


Since 

00 

r  /  -»  \  ^"*1  -i 

i  a-p)  p  =  i 

n=l 

the .probability  of  eventual  acceptance  is  one.  On  the  other  hand, 
the  expected  waiting  time  is 

CO 

„  z,  Nn-1  1 

E  n(l-p)  p  =■  - 

n=l  y 

From. this  it  is  clear  that  the  speed  of  Method  I3>  is  improved  by 
increasing  p.  This  can  be  done  by  first  of  all  choosing 

c  =  ^ 

'max 

Secondly,  (b-a)  can  be  decreased  by  truncating  any  long  thin  tails 
of  <P(x). 


Under  any  circumstances  Method  II  has  the  serious  disad¬ 
vantage  that  no  useful  upper  bound -on  its  running  time  can  be  given. 


2.9-b-  Generating- the  -Normal  Distribution 

Fortunately  Method  I  and  Method  II  do  not  exhaust  the 
possibilities  for  designing  random  number  generators.  As  -an  example 
of  the  exception  we  shall  now  describe  a  unique  method  for  generating 
random -numbers  from 'the  normal  distribution  with  mean  0  and  standard 
deviation-1.  In  this  case  Method  I  is  not  easily  applicable  because 
there  is.no  simple  formula  for  the  function  $(x),  and  Method  II  is 
undesirable  for  the  reasons  we  have  already  mentioned. 

First  of  all,  let  us  note  that  Method  I  yields  an  efficient 
random  number  generator  for  the  distribution  with  density  function 


Here 


and 


^_1(y)  =  (-lnCl-y))1^2. 


Furthermore,  if  we  choose  r  from  the  above  distribution  and  0  from  the 
uniform  distribution  ( 0 ,2tr ) ,  the  distribution  density  function  for  the 
pair  (r,0)  is  given  by 
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2ir 

No  if  suppose  that  (r,0)  are  the  polar  coordinates  of  a  point  in  the 
plane.  Then  if  we  change  to  Cartesian  coordinates 

x  =  r  cos  0,  y  =  r  sin  0 

the  density -function  becomes 

Because  the  density  factors  in  this  way,  we-  conclude  that  x  and  y 
are  independent  and  each  is  normally  distributed. 

2. 9.-5  Uniform  on -the  Unit  Sphere 

A  variant  of  Method  II  is  often  used  to  sample  from  the 
uniform  distribution  on  the  unit  sphere  using  the  following  procedure 

A.  Choose  three  independent  random  numbers  rj,  r2,  r3  from  the 
uniform  distribution  on  (0,1). 

B.  Retain-  if  0  <  Z  r.2  <  1;  otherwise  go  to  A. 

C.  The  required  unit  vector  is  (ri2+r22+r32)  fri *^2 »f3) • 


2.9 '.6-  Effusive  Flow 

< 

Finally,  we  describe  a  modification  of  Method  II  for 
generating  random  numbers  from  a  distribution  with  a  density  function 
proportional  to 

x  e-(X-Xo^y^:2  on  (0,b), 

Ideally  we  would  like  to  consider  shis  function  on  (0,»),  but  the 
cut-off  at  b-must  be  introduced  to  control  the  expected  waiting  time. 
This  particular  distribution  is  important  because  it  occurs  frequently  - 
in  dealing  with  boundary  conditions  for  the  Boltzman  equation. 

The  algorithm  is  as  follows: 

A.  Choose  r  from  the  normal  distribution  with  mean  zero  and  standard 
deviation  one .  Then 

x  =  x  +  r  c/21/2 
o 

has  density  function 


1  -(x-x  )2/c2 

-  e  p 

Cv/tT 
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B:  If  x.^  (o,b)  we  reject  this  value  of  x  and  repeat  step  1,. 
Otherwise* -we'proeeed  to  C. 

•  G;  Choose  a.  random  number  y  from  the- -uniform- distribution  on  (0,b). 
If  x  <  -yrwe  -accept  x  as  the  required  random. number.  Otherwise, 
the.  program-'goes.  to.  A  -to. .repeat-  the;  process . 

■We  shall  omit  the  details-of  the  -theory  in-  this-  case  since 
they  follow; -closely  our  discussion -of -Method 'll.  The -algorithm  has 
an'.unreasonably -large-expected -waiting  time- if  x  <  -c.  The  waiting 
time,  can- be -reduced  if -x  <  0  by  using’-only  the ‘positive  part  of  the 
.normal  distribution  -in  -s?ep  1,  i.e. 

x  =  x^+  |r | c/21^2 

If  x  is.  large- and--negative- then  the  standard  Method 'II  is  more 
efficient -than- this  one. 
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.3*0  -Introduction 

J.  K.  Kaviland  in  1961  introduced  a  Monte  Carlo  method  for 
obtaining- time  independent  solutions  of  the  -Boltzman  equation  which 
is .now  ^descriptively  called  the  “test  particle  method" -  This  method 
requires. enough. storage  for  two  approximations  to  the  distribution 
functions  .-to  .be  in  memory  at  -any  one  time  and  has  therefore  been 
-limited  to  problems  which  could  be  reduced  to  one  space  and  two 
velocity  dimensions:  The  computing  times  are  of  the  order  of  a  half 
hour  on  the  T.B.M.*  --709  to  -obtain  -a  density  profile  with  a  standard 
deviation  of. the  ’order  -of  1/20  -of  the  -density  at  any  given  point. 

In  this  case -the  density-histogram  was  accumulated  in  10  cells  each 
having,  a -width- -of ’the  order  of  1/10  mean  free  path. 

3.1  Preliminary  Theory 
3.1.1  Outline 


The  basic  idea  in  Haviland's  method  is  as  follows:  In 
order  to  compute  the  distribution  function  f(x,v)  satisfying  a 
given  Boltzman  equation,  and/or  some  of  its  moments  in  a  given 
3-dimensional  region  3)  with  given  boundary  conditions  on  the  boundary 
d&  of  3) .  •  Denote  the  phase  space  for  this  problem  by  ft  =  E3. 

Since  the  computation  is  designed  to  only  produce  a  histogram  approx¬ 
imating  the  distribution  function,  divide  the  phase  space  into  a 
finite  number  of  cells  ft^  (in  practice  k  will  represent  a  vector  of 
3  integers). 

Initially  one  makes  an  educated  guess  (usually  an  appropriate 
linear  combination  of  Maxwellians)  for  an  approximate  distribution 
function.  This  is  considered  the  target  distribution  for  the  first 
iteration.  The  target  distribution  is  then  used  to  statistically 
compute  trajectories  in  ft  of  a  molecule  having  the  same  collision 
cross-section  with  the  target  molecules  as  appears  in  the  given 
Boltzman  equation.  As  the  trajectories  are  computed,  their  contrib¬ 
ution  to  an  incident  distribution  function  is  recorded.  Clearly,  if 
the  initial  guess  was  an  exact  solution  to  Boltzman 's  equation,  one 
should  find  that  as  the  number  of  trajectories  increased  the  partially 
accumulated  incident  distribution  function  would  tend  to  the  target 
distribution  function.  Otherwise,  the  incident  distribution  approx¬ 
imation  is  different  from  the  target  distribution  even  after  a  large 
enough  number  of  the  trajectories  have  been  computed  to  reduce  the 
variance  to  an  acceptable  level. 

Nevertheless,  it  has  been  found  by  experience  that  the 
incident  distribution  is  a  better  approximation  to  the  correct  answer 
in  the  sense  that  if  the  process  is  repeated  with  the  target  distribution 
replaced  by  the  previously  computed  incident  distribution  then  eventually 
the  incident  distribution  will  tend  to  the  target  distribution.  This  is 
interpreted  as  convergence  of  the  iteration  scheme  to  the  correct  dist¬ 
ribution  (=  final  incident  distribution  =  final  target  distribution). 

The  required  number  of  iterations  depends  on  the  initial  guess  and  the 
error  tolerances.  Haviland  used  fewer  than  5  in  his  calculations. 
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3-1.2  Trajectories  to  Distribution  Function 

To  earolain  and  motivate  Havi land's  procedure  for  accumulating 
the  incident  distribution  function,  suppose  that  «&,  ft,  ft^,  have  the 
same  meaning  as  before.  Imagine  that  the  molecules  in  «D  are  numbered 
with  an  index  a  (a  =  1,2,..., H),  and  have  trajectories  given  by 


and 

-  - 

x  (t)  =  x  (0)  +  I  v  (t)  dr. 

a  a  J  a 

J  o 

At  this  point  it  is  worthwhile  to  remark  that  each  v  (t)  must  closely 
approximate  a-  stepfunction  (piecewise  constant  function)  in  order 
that  the  Boltzman  eauation  be  valid  for  this  system. 

The  first  task  is  to  define  what  is  meant  by  the  statement 
"f  is  the  distribution  of  the  dynamical  system  whose  trajectories  are 
x  (t),  a  =  1,2,. ..,N": 


6  (x  (t),v  (t))  +  R 
k  a  a 


(1) 


where  6^  is  a  function  having  the  property  that  Sj^Xjv)  =  1  when 
(x,v)  belongs  to  the  cell  ftk  and  5k(x,v)  =  0  otherwise;  and  R  is  a 
random  variable  having  mean  zero  and  variance  tending  to  zero  as 
the  volume  of  and  therefore  N  becomes  large.  Strictly  speaking, 
one  should  define  the  probability  space  underlying  R  but  it  is  best 
to  avoid  such  a  project  at  this  point  for  the  sake  of  brevity.  In 
fact,  in  order  to  elucidate  the  essentials  R  shall  be  omitted 
henceforth  and  =  replaced  with  =  to  indicate  the  omission. 


Since  (l)  is  true  for  all  t  it  follows  that 


T  s'  NT 

I  /  Ja  ^  I  (2> 


It  is  easy  to  see,  however,  that 

-T 


f  Sk(xa(t),va(t))dt 

w  r\ 


is  nothing  other  than  the  transit  time  t  ,k  that  the  trajectory  of 
the  a-th  molecule  spends  in  the  k-th  celS.  The  factor  —  is  the  same 
for  each  cell,  and  it  is  more  efficient  in  actual  computation  to 
simply  accumulate  the  fc-s  in  a  suitable  array  and  then  divide  by 
the  normalizing  factor  at  the  end  of  the  computation.  Finally,  there 
is  a  very  important  observation  to  be  made  at  this  point.  It  is  not 
really  necessary  to  compute  the  actual  trajectories  of  N  molecules 
in  order  to  apply  the  preceding  formula.  It  suffices  to  generate 
trajectories  which  have  the  same  probability  distribution  as  the 
actual  ones.  It  must  be  noted  that  to  precisely  specify  this 
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distribution  would  entail  a  certain  amount  of  basic  work  which  does  not 
yet  appear  in  -the  literature .  One  must  he  satisfied  with  describing  a 
generation  procedure  which  nicies  statistically  the  one  occurring  in 
nature  and  arguing  by  analogy  rather  than  proof. 


3.1-3  Computation  of  Trajectories 

Suppose  a  molecule  starts  at  (x  ,v  )  in  phase  space  at  time 
t  =  0.  Then  the  probability  that  the  molecule  will  undergo  its  first 
collision  in  the  time  interval  (t,t+dt)  is  given  by 

P(xq,  vq,  t)  =  (Z  e-tS)dt 

where 

C  ^b  ,2.  u 

Z(x  ,v  )  =  Idv  I  m  bdb  j  de  g  f  (x  ,  v)  (3) 

0  0  j  J  O  J  o  0 

is  the  collision  rate  with  the  target  molecules.  As  is  customary, 

g  =  |v  I  ,  and  v  =  v  -  v  . 
r1  r  o 

Let  t  be  chosen  at  random  from  a  distribution_with  density 
P(x  ,v  ,t).  The  position  at  time  of  collision  is  then  xj  =  x  +  tv  . 
In  8rder  to  obtain  a  post  collision  velocity  for  the  molecule°reason 
as  follows:  The  collision  partner  has  a  velocity  chosen  from  a  dist¬ 
ribution  with  density 


"pfxT  =  dv )  (*0 

The  probability  that  the  collision  parameters  b,e  fall  into  the  intervals 
(b,b  +  db),  (e,e  +  de)  is  given  by 

bdbde  ^ 

irb2 

m 

Let  them  be  chosen  accordingly.  Then  the  post  collision  velocity  V 
may  be  computed  using  (see  1.1.3) 

V =  v +  ( v  • e )  e  (6) 

00  r 

3.1.1*  Specialization  to  One  Space  and  Two  Velocity  Dimensions 

Haviland  considered  two  problems  which  require  only  one  space 
dimension  and  two  velocity  dimensions  for  representing  the  distribution 
function:  the  flow  between  infinite  parallel  walls  at  two  different 
temperatures  and  a  plane  shock  extending  to  infinity  in  both  directions. 
In  both  these  problems  the  distribution  function  can  be  immediately 
written  f(x,v  ,v  )  (see  1.2.5). 

X  "*■ 

Let  the  (x,u,v)-space  be  partitioned  into  cells  having  dimen¬ 
sions 

Ax,  Av  ,  Ay;. 
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The  number  of  molecules  .in  a  typical  cell  is 


/x+nx  sv  +Av-  rXi.'ra'tL 

dx"-J  X'  Xdv' J  &v±  J  T^d0f(x,vx,v, 


V,_+AVj_  y-2v 


=  2tt  f (xjV^,-^)  v^AxAv^Av^ 

molecules  per  unit  area  where  (x,v  •  3Vj_ )  is  some  interior  point  of  the 
cell.  We  now  define 

FUjV^VjO  =  2n  vx  f(x,vx,vx).  f7) 

Suppose  each  cell  in  phase  space  is  now  indexed  with  subscripts  (i,j,k) 
giving  the  coordinate  numbers  of  the  cell  in  the  x,vxand  directions 
respectively,  and  (x.  ,v  ,vj^)  gives  the  midpoint  coordinates  of  the 
i  j  k-th  cell.  If *F^  gives  the  accumulated  trajectory  residence  time 
in  the  ijk-th  cell,. Men 

F(x.  ,v  .  ,v  -  )  =  -  F...  (8) 

1  xj  J-K  AxAv  Av,  ljk 

x  J-  u 

where  C  is  a  normalization  factor  chosen  so'that  the  distribution  function 
predicts  the  correct  average  number  density  of -molecules.  In  order  to 
compute  collision  rates  and  moments  of  the  distribution,  the  following 
formula  given  by  Haviland  is  very  useful 

|f(x,v  ,v  ,v  )  $  (x,v  )dv  dv  dv 
J  *  x’  y*  z'  *  x  x  y  z 

/“  ^2tt 

dv.j  j  v|_d0f(x,v  J-^cosOjV^sinO)  4>  (x,v  ,v^cos0,v^ sin©) 

_«  XJ  o  J  o  X 

,  s-™  ^27t 

=  j  dvx|  dv  F(x,v  jVj^)  j  d©  $  (xjV^^cosOjV^sin©)  (9) 

i  r2lT 

!  ar  A?i»  J  a0  4  <xi-\j’v-Lkcose’vxksin£,) 

J  ,K  O 

where  x^  is  the  closest  cell  midpoint  to  x. 


3.1.5  A  Finite  Velocity  Range  is  Sufficient 

Since  it  is  expected  that  for  large  values  of  v  ,vx  the 
distribution  function  will  decay  like  exp-[(vx/c  )2+(vL/cx)2] ,  only 
the  cells  near  the  mean  velocity  should  contain  ,  any  significant  num>  „-r 
of  molecules.  This  leads  to  the  very  important  simplification  of  working 
with  only  a  finite  number  of  cells  in  velocity  space.  Haviland  (p.151, 
table  IV)  has  studied  the  effects  of  restricting  the  velocity  distribution 
to  a  rectangle  of  the  form 

(u  -ac  ,u  +ac  )  X  (u -ac  ,u.+ac  ) 

X  o7  X  O  X  o  -*•  o 
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(where  (ux,ux)  is  the  mean  velocity  and.c0  is  the  most  prohahle  deviation 
from- the  mean)  for  various  values  of  a.  He  points  out  that  if  a  =  2.0 
then-  one-  is  assured  of  including  over  99%  of  the-  molecules  of  an  equil¬ 
ibrium  distribution.  Any  molecules  found  outside  the  rectangle  are  simply 
assigned  to  the  closest  cell  on  the  boundary  of  the  rectangle. 


3. -2  The  Algorithm 

3.2.1  Storage .Requirements 

The  largest  arrays  are  the  target  distribution  function  FT; 
the  incident  distribution  function  FI  each  of  dimensions  I,J,K,  a  cross- 
section  .array  XSECT  of  dimensions  K,  (.2J  -  l),  K  and  a  collision  rate 
array  SIGMA  of  dimensions  I,J,K.  Both  the  cross-section  and  collision 
rate  arrays  are  used  to  save  time  by  storing  functions  which  would 
otherwise  have  to  be  recomputed  a  large  number, of  times.  Haviland  used 
I  =  10,  J  =  2U,  K  =  ,12  in  all  of  his  calculations. 


3.2.2  Data  and  Formulas 


In  the  case  of  flow  between  parallel  plates  at  different 

temperatures  with  perfect  accommodation  at  the  walls,  there  are  two 

non-dimensional  parameters  -  the  temperature'  ratio  of  the  walls  and 

the  distance  between  them  in  mean  free  paths.  This  is  also  the  stage  at 

which  any  constants  which  will  be  used  frequently  during  the  program 

should  be  r'ead  in  or  computed.  They  include  the  number  rmax  of  iterations 

necessary  to  get  convergence  of  the  distribution  function,  the  number 

nmax  of  cell  transits  necessary  to  accumulate  the  next  approximation 

to  the  distribution  function,  and  the  array  XSECT.  In  Haviland' s 

calculations  r  '  3  and  n  ~  50,000-100,000. 
max  max 

The  formulas  for  both  XSECT  and  SIGMA  are  derived  from 


/oo  co  2lT 

dv'  f  dv£  F(x,v',vj0  f  d©gb2 

_ao  r\  %j  r\ 


o  -  o  m  (10) 

which  is  obtained  by  using  (3),  (T)s  and. (9).  The  evaluation  of 

■  2ir 

)cl 

m 


J  de< 


at  the  cell  midpoint  values  (v^,^^  ana  (vxj»XLk^  y*elds 

XSECT  (k,J'-J+J,k') 

=  1/2  f  b2  [(v  .-v  /)2+(v;\  cosO'-v. ,  sin0)2+(v/\  sin0'-v  .  sin0)2]^^ 
J  m  xj  xj  Ik  lk  Ik  xk 


since 


A  couple  of  remarks  pertain  to  the  preceding  formula.  Firstly, 


(vx^os0r-vXkCos0)2+(v^sin0'-vikSin0)2^i;2+vL2-2v;kV'J_kcos(0'-0)  ) 
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and  the  integration  is  over  a. full  period  of  the  cpsine  function,  the 
integral  is  independent  of  O'.  Secondly,  since  the  jells  in  the  vx 
direction  are. equally  spaced,  vxJ-vXj  will  depend  only  on  j'-j.  This 
is  the  justification  for  writing 'XSECT  (k,j‘'-j+J,k")  instead  of  XSECT 
(k,j  ,j',k'):  The  J  has  "been  thrown. in  for  the  .sole .purpose  of  keeping 
indices 'positive ■ -  a  requirement  of  the  FORTRAN  language  in  which  this 
program  was  originally  written.  It  is  also  worthwhile  to  note  that 
SIGMA  should  he  computed -after  FT  which  changes  with  each  iterate. 

Initially,  the  starting  target  distribution  function  is  stored 
in  FT.  'For  this  purpose,  -it  .is  a  good  idea  to  use  the  best  analytical 
approximation  which  is  known  and  reasonably  simple.  Haviland  used  the 
free  molecule,  solution  for  the  heat  transfer  problem  and  the  bimodal 
Mott-Smith  distribution  in  the  shock  problem.  If 


f(x,vx,vjcos0,v^sin0)  =  fCxjV^v^ 
is  the  initial  analytic  approximation,  we  begin  by  setting 


[Vta/2  fVj.^/2 

ij,k)  =  2tt  dx  dv  I  dvavxf(x,v  ,vA) 

x, -Ax/2  J  v  ,-Av  /2  Jyt,-AyL/2 


FT(i 


"v  ,+Av  /2 
■  x 


vik+AvJ./2 


x.^-Ax/2 


,  .  -Av  /2 
xj  x' 


Y1_k-Av-a./2 


for  the  interior  cells.  For  boundary  cells,  the  tail  of  the  distribution 
must  be  included.  For  example,  if  j  *  J  above  then 


is  replaced  by 


hv  +Av  /2 
x  x 

v  v  -Av  /2 

X  X 


/CO 

VV2 


3.2.3  Normalization 


The  normalization  constant  Cp  is  computed  from  the  formula 

n  Ax  =  C  Z  F(i,j,k) 

In  r  . 

ijk 

where  is  the  mean  number  density.  Note  that  the  left  side  of  this 
formula  is  nothing  other  than  the  total  number  of  molecules  per  unit 
area  (since  the  system  is  supposed  to  extend  tc  »  in  the  y  and  z 
directions).  Any  moments  which  one  wishes  to  compute  and  output  can  be 
dealt  with  simply  once  the  normalization  constant  for  the  distribution 
function  is  known. 


3.2.4  Trajectory  Initialization 

The  first  task  which  must  be  performed  even  before  any 
trajectories  are  considered  is  the  computation  of  the  array  SIGMA 
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according  to  the  formula 


SIGMA  (i,j,k)  =  (C  /Ax)  Z  FT(i,j',K')  XSECT  (kJ'-j+J,^)  • 

F  r,k' 

The  actual  computation  of  a  trajectory  starts  with  the  molecule 
at  either  the  left  boundary  or  right  boundary  •  having  velocity  (vx,vx) 
c-hosen  according  to  the  boundary  conditions  and  is  repeated  every  time 
the  molecule  crosses  a  boundary.  For-  many  kinds  of  boundary  condition 
the  vx,vj_  are  chosen  from  an  appropriate  probability  distribution.  The 
matter  of  boundary  conditions  is  a  study  in  itself,  and  the  present 
discussion  is  restricted  to  equilibrium  .flow  normal  to  a  plane  and  to 
emission  from  a  perfectly  accommodating  wall. 

Let 

$(u,c;vx,vx)  =  C$vx  exp(  -  ((vx-u)2  +  vf  )/c2) 

where  C$  is  chosen  so  that  /$  dvx  dvx=  1,  and  c  is  the  most  probable 
speed  after  the  mean  is  subtracted  out.  In  the  case  of  a  flow  in  the 
x  direction  of  speed  u,  the  probability  that  a  molecule  which  crosses  a 
fixed  plane  normal  to  the  flow  has  velocity  in  the  range  vx,vx+dvx; 
vx  ,VjL+dvL  is  given  by  $(u,c;  vx,vx)  dvxd\Q_ .  In  the  case  when  the  molecule 
is  emitted  from  a  perfectly  accommodating  wall  the  probability  is 
4>(0,c ;vx,Vj_ )  dvx  dvj_.  Efficient  procedures  for  choosing  vx,v  at  random 
from  a  distribution  with  density  function  $  have  already  been  described 
(2.9.6). 


3.2.5  Cell  Coordinates 

Once  the  velocity  components  of  the  test  molecule  are  chosen, 
its  cell  indices  i , j ,k,  are  computed  for  future  reference,  but  the  actual 
coordinates  are  not  changed  to  the  cell  midpoint  values.  The  molecule 
now  has  two  choices:  to  cross  a  cell  division,  (crossing  the  calculation 
field  boundary  is  a  subcase  of  this)  or  to  undergo  a  collision  within  its 
cell.  In  order  to  decide  which  occurs,  the  program  computes  tc,  the  time 
to  next  collision  and  the  time  to  cross  a  cell  division.  The  expression 
for  follows  from  a  simple  geometrical  exercise. 

As  a  basis  for  discussing  the  computation  of  tc  refer  to  (lQ). 
The  array  SIGMA  has  already  been  computed  and  the  indices  i , j ,k  are 
computed  in  this  section.  Choose  tc  from  the  distribution  whose  density 
is 

SIGMA  (i,j,k)  e-t  SIGMA 


using  Method  I. 

If  tc  >_ 2.6  is  applied  to  move  the  molecule  into  the  next 
cell.  Otherwise,  2.7  is  used  to  compute  a  collision. 


3.2,6  Translation 

The  x  coordinate  of  the  molecule  is  changed  to  the  appropriate 
cell  division  value  and  is  added  to  FI  (i,j,k).  If  the  molecule 


Ul 


happened  to. end  up  at  the  edge  of, the  x  interval  2.5  is  applied  to  start 
a  .new  trajectory.  Otherwise,  the  collision/transit  calculation  is -repeated 
in  the  new  cell. 


3.2.7  Collision 

To  start  with,  t.->  is  added  to  Fl(i,j,k).-  Then  the  cell 
(indices  j"jk')  where- the  collision  partner  will  originate  is  chosen  in 
accordance  with  the  unnormalized  probability  function 

FT(i,j',k")  XSECT  (k, j '-j+J,k') . 

Although  Haviland  used  Method  I  to  make  his  choice  of  j",k''  it  seems 
that  in  this  case  Method  II  would  be  faster,  even  though  it  would  entail 
finding  the  maximum  values  of  Ff  and  XSECT  in  2.4  in  order  to  determine 
a  suitable  range  for  the  second  random  number. 

Once  the  choice  of  yk'  is  made,  the  collision  partner  is 
assumed  to  have  velocity  components  (v^,v£)  equal  to  the  cell  midpoint 
velocity  coordinates.  In  order  to  obtain  cartesian  components,  choose 
0'  from  the  uniform  distribution  on  (0,2tt)  and  let 

v'  -  v'cosQ v'  =  vf  sin©'  . 

y  X  »  z  x 

The  most  efficient  way  to  compute  the  post  collision  velocity 
V  of  the  incident  molecule  depends  to  a  large  extent  on  the  intermolecular 
interaction.  In  the  case  of  hard  spheres,  where  it  is  known  that  the 
post-collision  relative  velocity  (1.1.3)  has  direction  uniformly  dist¬ 
ributed  on  the  surface  of  the  unit  sphere,  it  is  fastest  to  choose  a 
unit  vector  ft  according  to  2.9.5  whereupon 

V  =  |v'  — v I w 
r  1  1 

and  (see  1.1. l) 


V  =  l/2(v  +  v'-  V  ) 

In  general,  when  the  distribution  of  w  is  not  expl-icitely  known, 
the  general  theory  of  Chapter  1  must  be  used  to  compute  collisions. 


3.2.8  Exit  Criteria 

The  collision-transit  computation  (2.5,  2.6,  2.7)  is  repeated 
until  the  total  number  of  collisions  plus  transits  reach  a  preassigned 
number  sufficiently  large  for  the  required  standard  deviation  in  the 
results  (see  3.3)* The  attainment  of  this  number  signals  the  end  of  the 
iteration  and  the  beginning  of  the  next  one  or  the  completion  of  the 
calculation . 

In  proceeding  to  the  next  iteration  FT  is  stored  in  FI.  The 
need  for  an  actual  storing  operation  can  be  eliminated  by  a  suitable 
arrangement  of  the  program  and  it  is  test  to  refer  to  Haviland' s  original 
flow  diagrams  to  see  how  this  is  done.'  ' 


This  completes  the  description  of  the  algorithm. 


3.3  Variances 

3.3.1  Variance  of  FI  (i,j,k)  •< 

Haviland's  method  lends  itself  easily  to  a  theoretical  estimate 
of  the  variance  of  the  answers  it  produces.  .  In  order  to  obtain  this 
estimate,  Haviland  did  not  analyse  his  algorithm  directly  hut  instead 
assumed  a  model  which  is  supported  by  the  conclusions  it  produces  when 
the  distribution  function  is  Maxwellian.  The  model  is  along  the  .following 
lines:  Suppose  that  each  FI,(i,j,k)  is  the  sum  of  N  independent  random 
variables  Xa  each  taking  the  value  t  with  probability  F^j^  and  0  with 
probability  1  -  where 


Then 


and 


1  fnk *  1- 
ijk  i,)k 

E(V  “  TFiJk 

Var  (X  )  =  E  <X|)  -  W*,)')2  =  T2  F.Jt  (1  -  F.^) 


It  follows  that 


and 

(See  2.2.7) 

In  practice 


E(F(i,j,k))  =  E(EXa)  =  NtF 
a  J 

Var  (EXa)  =  Nt2  F±^  (l  -  F^). 


F.  ..  =  Fl(i,j,k)/  E  Fl(-i,j,k) 
ijk 

is  used  as  an  estimate  of  Fijfc.  For  reasons  of  simplicity  however, 
suppose  t  is  known  and  one  can  therefore  use 

Then 


as  required  and 


E(F. ,,  )  =  F. 
ijk  ijk 


<FiJk>  *  ?  Fijk  (1  -  Fi3k> 


3.3.2  Variance  of  the  Moments 

If  <Mx,vx,Vj. )(or  ^(vx,v^)  is  a  given  function,  then  the  moment 
of  ip  with  respect  to  the  distribution  F(x,vx,v,i,)  is  defined  to  be 


b3 


/^(x5vx,v j.)  ?{x.vy,vj)  ax  avx  av^ 

(or  f'H  vx,V|_)  ?(x,vx,Vj.)  dvx  qv_l  respectively) 

Consider  only  the  first  case  which,  in  its  discrete  version  * 
will  turn  out  to  contain  the  second-  The  discrete  approximation  for  the 
moment  of  $  is  given  by 


#  =  2 


p 

Vi,j,k.  ijk 


and  is  again  a  random  variable.  Then 


i,J.k 


'^ijk  Fijk 


Var  (ijO  =  ^ 


ti)2  p  (T_p  ) 

vijk  ijk  u  ijk; 


This  is  not  the  same  formula  for  the  variance  of  as  derived 
by  Haviland,  because  the  initial  model  was  slightly  different.  It  has 
the  advantage  of  predicting  a  larger  variance  than  Haviland' s  formula; 
this  seems  to  be  desirable  in  view  of  the  fact  that  Haviland 's  theoretical 
variances  turned  out  to  be  smaller  than  his  experimental  (computed)  ones. 


3.1*  Conclusion 

Haviland 's  method  was  the  first  attempt  to  solve  the  Boltzman 
equation  by  a  Monte  Carlo  technique  and  consequently  has  played  an 
important  role  in  the  development  of  the  subject.  Its  storage  and 
computing  time  requirements  pretty  well  limit  its  applicability  to  the 
two  types  of  problem  originally  treated  by  Haviland.  It  does  have  the 
advantage,  however,  of  being  amenable  to  a  close  mathematical  analysis 
(which  at  the  time  of  writing  has  not  been  fully  carried  out)  yielding 
a  proof  that  statistically  computed  quantities  have  the  required  mean, 
and  estimates  of  the  variances  of  these  quantities.  It  is  fairly 
reliable  as  such  methods  go,  and  so  far  the  only  unexplained  predictions 
have  been  regular  density  fluctuations  in  the  high  temperature  high 
density  regions  of  a  flowfiela. 
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4.  M0KTE  CARLO  SIMULATION  METHOD 


4.1  General  Description 

The  norsnalizability  (1.3.1)  of  the  Boltzman  equation  suggests 
that  an  actual  flow  involving  in  the  vicinity  of  1020  molecules  could 
he  studied  by  means  of  a  model  using  a  small  number  of  molecules,  each 
of  proportionally  larger  cross  section.  There  are  a  number  of  diffic¬ 
ulties  which  arise  from  such  a  trade-off.  The  higher  probability  of 
multiple  collisions  must  be  dealt  with  in  the  model  and  the  statistical 
fluctuations  arising  from  the  small  sample  size  necessitate  averaging 
a  large  number  of  observations  to  obtain  a  useful  approximation  to  the 
distribution  function. 

Nevertheless,  G.A.  Bird  and  others  have  had  considerable 
success  using  a  model  with  as  few  as  1000  molecules  in  conjunction 
with  .a  Monte  Carlo  procedure  for  obtaining  post  collision  velocities 
consistent  with  the  cross-section  function. 

Bird's  method  is  an  algorithm  in  which  .each  iteration 
operates  on  a  sample  of  N  (of  the  order  of  a  thousand)  molecules 
from  a  distribution  with  distribution  function 

f(x,v,t) 

to  generate  a  sample  from  a  distribution  with  density  function 

f(x,Vjt+At) 

‘where  f  satisfies  the  Boltzman  equation  and  suitable  boundary  conditions. 
Although  the  foregoing  statement  has  never  been  precisely  proved,  the 
method  has  found  extensive  application  in  the  study  of  problems  in 
the  transition  regime  wher.  the  mean  free  path  is  of  the  same  order  of 
magnitude  as  the  flow  properties  under  study.  In  almost  all  cases,  the 
results  are  quite  reasonable.  The  method  is  applicable  in  both  time 
dependent  and  time  independent  problems.  In  the  case  of  time  dependent 
problems  it  is  necessary  to  repeat  the  calculation  using  different 
random  number  sequences  a  sufficient  number  of  times  to  reduce  the 
standard  deviation  of  any  sampled  quantities  to  an  acceptable  level. 

In  the  case  of  time  independent  problems  the  usual  procedure  is  to 
start  with  an  arbitrary  distribution  function  f(x,v,o)  and  to  run  the 
iteration  without  sampling  until  f  has  converged  sufficiently  to  its 
steady  state  value  f(x,v,“).  Then  sampling  is  done  at  regular  inter¬ 
vals  until  a  sufficient  sample  size  is  reached. 


4.2  The  Algorithm 

4.2.1  Storage  Layout  and  Initialization 

Each  molecule  is  represented  in  computer  memory  by  storing 
its  velocity  and  position* coordinates .  This  requires  3+3=6 
locations  but  in  many  cases  geometrical,  symmetry  allows  a  reduction 
in  the  number  of  position  coordinates  which  must  be  stored.  Many 
physically  interesting  problems  (for  example,  the  flow  around  spher¬ 
ically  of  cylindrically  symmetrical  obstacles)  require  only  two 
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position  coordinates;  and  some  studies  on  the  relaxation  of  the  distribution 
require  no  position  coordinates.  Space  is  assigned  in  the  fast  access 
memory  of  the  computer  to  accommodate  the  coordinates  of  up  to  simulator 
molecules . 


There  remains  the  problem  of  referencing  the  molecules  in  such 
a  way  as  to  make  those  in  an  arbitrary  region  of  position  space  easily 
accessible.  This  is  accomplished  by  dividing  position  space  into  cells. 
Some  notation  will  be  helpful  later.  Let  the  flow  region  under  study  be 
denoted  by  <>D  and  subdivided  into  a  finite  number  K  of  subregions  <2^ 
which  will  be  called  cells.  The  cells  must  be  small  enough  so  that  the 
velocity  distribution  function  is  approximately  constant  on  each  one; 

This  judgement  must  be  made  before  the  solution  to  the  problem  is  known, 
but  experience  has  shown  that  intuition  is  an  adequate  guide  in  setting 
up  the  cell  size.  Factors  affecting  the  choice  of  cell  size  are  considered 
in  Section  3.3. 


Once  the  cell  boundaries  have  been  decided  the  programming 
problem  of  easily  referencing  the  subset  of  molecules  in  each  cell  can 
be  tackled.  There  are  several  well  known  ways  of  doing  this  (refer  to 
any  work  on  list  processing)  of  which  one  possibility  is  to  set  up  K 
auxiliary  tables  ,in  such  a  way  that  the  k-th  table  contains  the  add¬ 
resses  or  sequence  numbers  of  the  molecules  in  the  k-th. cell.  The  method 
described  presently  combines .access  speed  with  storage  economy. 

Suppose  cell  k  contains  %  molecules,  1  <_  k  <_  K.  Then  the 
molecules  from  cell  1  are  stored  in  locations  1  to  ni  of  the  molecule 
listi  from  cell  2  in  locations  n\  +  1  to  n2  etc.  A  vector  MAP  with  K 

entries  containing  the  numbers  nj ,  nj  +  .  ni  +  n2  +  ...  +  n^ 

is  also  stored.  Whenever  an  addition,  deletion  or  move  occurs,  only 
the  last  molecule  in  each  cell  section  is  moved  and  MAP  is  updated  to 
describe  the  new  configuration. 

Preliminary  to  starting  the  iteration,  a  number  N0  <  Nmax, 
of  molecules  is  chosen  from  an  appropriate  starting  distribution  and 
stored  in. the  molecule  table,  of  the  programme.  Any  additional  information 
necessary  to  reference  the  molecules,  a  table  of  cell  volumes,  cross- 
section  constants,  boundary  condition  parameters,  etc.  must  also  be 
assembled  at  this  stage,  The  appropriate  starting  distribution  depends 
on  whether  the  algorithm  will  be  applied  to  a  steady  pr  unsteady  flow 
prob3ciu  and  is  discussed  in  Section  3.1. 


*+•2.2  Steps  -Comprising  One  Iteration 

A.  Collision  Inqrement  to  the  Distribution  Function 


The  following  procedure  is  repeated  to  each  cell  k,  1  £  k  <_  K. 
A  pair  of  molecules  is  chosen  at  random  from  those  in  cell  k  in  such  a 
way  that  the  probability  of  choosing  a  given  pair  is  proportional  to  the 
relative  velocity.  To  avoid  computing  l/2nk  (nk-l)  relative  velocities, 
it  is  convenient  to  estimate  the  largest  relative  velocity  smax  which 
can  occur  with  any  frequency  and  to  make  the  choice  as  follows:  Choose 
two  independent  integers  from  the  uniform  distribution  on 
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k-l  k. 

(l  +  I  n. ,  E  n. ), 
1  1  1 


rejecting  the  pair  and  repeating  the  choice  if  the  two  integers  turn 
out  to  he  equal,  and.  compute  the  relative  speed  |v  |  =  g  of  the 
accepted  pair.  Choose  a  random  number  r  from  the  uniform  distribution 
on  (o,smax). ,  accept  the  pair  for  collision  if  r  <  g  and  otherwise 
reject,  repeating  the  procedure  from  the  choice  of  two  random  integers. 
As  soon  as  a  pair  is  accepted  for  collision  the  collision  aprameters  are. 
chosen  at. random  and  the  velocity  vectors  of  the  pair  in  question  are 
replaced  with  post  collision  velocity  vectors  computed  on  the  basis  of 
the  collision  parameters.  Details  of  this  calculation  are  given  in 
Ch.  1.1  (see  also  Ch.  2.9.5). 


After  each  collision  a  time  counter  t  is  advanced  according 
to  the  formula 


a  T0lk 
A"t2* 


where  A  is  the  cross-section  per  simulator  molecple,  and  vol^  is  the 
volume  of  cell  k.  Collisions  are  computed  until  the  total  time  advance 
exceeds  a  preset  value  Atm.  Both  the  time  advance  formula  and  the 
value  of  Atm  raise  a  number  of  delicate  questions  and  are  discussed  in 
sections  3-5  and  3.3  respectively.  Once  the  collision  calculation  has 
been  completed,  the  programme  proceeds  to  B. 


B.  Convection  Increment  to'  the  Distribution  Function 

— . t . — — >  .  ... - — T t-  ■■  ■-  .  ■  ■■  - - 

Each  molecule  is  moved  to  a  new  location  corresponding  to 
the  time  interval  Atm.  This  means,  of  course,  that  its  position 
coordinate  x  is  replaced  with  x  +  vAtm  in  case  of  cartesian  coordinates. 
It  may  happen  that  this  displacement  of  the  molecule  takes  it  outside 
the  region  £2,  across  a  cell  boundary,  or  across  the  boundary  of  some 
obstacle  in  the  flow.  In  each  of  these  cases  an  appropriate  strategy 
must  be  applied.  For  instance,  crossing  an  exterior  boundary  entails 
deletion  of  the  molecule  from  the  list;  crossing  a  cell  boundary 
necessitates  appropriate  changes  in  the  bookkeeping  tables;  and  coll¬ 
ision  with  an  obstacle  entails  replacement  with  a  molecule  chosen  from 
a  distribution  appropriate  to  the  boundary  condition  (see  Section  3.6). 
For  example,  the  molecule  may  be  reflected,  emitted  diffusely,  or  it 
may  choose  one  of  these  fates  with  probability  a.  When  the  whole 
molecule  list  has  been  dealt  with,  the  programme  proceeds  to  C. 


C.  Influx  of  Molecules 


Facilities  must  be  provided  in  the  programme  for  introducing 
molecules  into  the  calculation  which  enter  the  region  because  of  the 
flow.  These  are  usually  introduced  with  their  position  coordinates 
consistent  with  entry  into  the  region  at  some  random  time  in  the 
interval  (t,  t  +  Atm) ,  and  velocity  coordinates  appropriate  to  the 
incoming  distribution.  The  boundary  point  of  entry  is  chosen  at  random 
when  no  better  criterion  is  available  and  the  number  of  molecules 
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introduced  in  this  way  is  derived  from  the, theoretical  input  flux  and 

the  elapsed  time  At  (details  in  Section  3.6). 

m 

D.  Data  Collection 


The  description  to  this  point  constitutes  one  iteration. 
References  to  data  collection  have  been  omitted  purposely  in  order  to 
make  the  description  of  the  main  steps  as  concise  as  possible. 
Basically,  data  collection  falls  into  two  categories;  estimation  of 
various  transfer  rates  across  specifies  boundaries  in  the  flow,  and 
estimation  of  moments  of  the  distribution  function  at  fixed  time. 
Sampling  for  the  first  category  is  most  easily  carried  out  as  part  of 
step  B,  and  for  the  second  category  at  this  point  in  the  program.  The 
latter  is  .not  normally  done  after  every  iteration  and  the  actual 
strategy  is  dependent  upon  whether. a  steady  or  unsteady  flow  is  being 
simulated  as  described  in  Section  3.1. 


E.  Exit  Criterion 


When  the  programme  has  reached  this  point ,  two  alternatives 
are  possible:  commence  another  iteration  at  A  or  prepare  and  output 
final  data.  The  decision  can  be  made  either  on  the  basis  of  the  total 
number  of  iterations  or  the  accumulated  time  advance.  The  number  of 
iterations  necessary  to  accumulate  the  required  data  depends  on  the 
acceptable  standard  deviation' of  the  answer  (see  Section  3.7). 


U.3  Analysis 

U.3.1  Steady  and  Time  Dependent  Problems 

The  treatment  of  steady  and  time  dependent  problems  differs 
in  the  choice  of  the  initial  distribution  function  and  in  the  sampling 
strategy.  In  the  case  of  steady  problems  a  good  approximation  to  the 
actual  distribution  function  is  seldom  known.  Consequently  the  initial 
distribution  function  is  seldom  known.  Consequently  the  initial  dist¬ 
ribution  for  the  iteration  is  usually  taken  to  be  a  drifting  Maxwellian 
of  approximately  the  right  velocity,  temperature,  and  density.  This 
initial  distribution  is  expected  to  converge  to  the  final  steady  state 
distribution  fairly  quickly.  An  early  study  done  by  Bird  [.3  )  suggests 
that  convergence  takes  three  to  four  collision  times.  Sampling  is 
delayed  until  convergence  has  occurred  and  then  carried  out  at  regular 
intervals  until  a  sufficient  number  of  observations  has  been  made. 

If  the  sampling  interval  is  too  short  then  the  successive  observations 
are  not  independent  and  consequently  little  is  gained  in  the  way  of 
reduced  standard  deviation.  On  the  other  hand,  if  the  sampling  inter¬ 
val  is  too  long,  there  is  a  waste  of  computer  time.  The  proper  choice 
of  interval  depends  on  many  factors  (see  ^.3.7)  but  generally  speaking, 
it  has  been  found  that  1/U  to  1  collision  time  is  a  satisfactory  range. 

In  the  case  of  time  dependent  problems,  the  initial  distribution  function 
is  assumed  known  and  can  therefore  be  used  as  a  source  distribution  for 
the  initial  sample  of  molecules.  Sampling  is  done  at  those  times  when 
a  snapshot  of  the  flowfield  is  desired,  and  the  iteration  is  continued 
until  the  flow  has  developed  sufficiently.  In  order  to  reduce  the 
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standard  deviation  of  the  results  to  an  acceptable  level  it  is  nec¬ 
essary  to  repeat  the  calculation  a  large  number  of  times  using  a 
different  random  number  sequence  each  time,  and  to  average  corres¬ 
ponding  results. 


h.3.2  Cell  Size  and  At 

m 


From  comments  in  the  literature  [,  6,  Fig.  2.]  and  from  the 
author's  own  experience  it  seems  that  a  vide  range  is  available  for 
both  the  cell  diameters  and  the  time  .advance  interval  without  signif¬ 
icantly  altering  the  results.  There  are  extremes,  however  which 
should  be  avoided.  Since  the  method  treats. the  distribution  function 
as  if  it  were  constant  over  each  cell,  the  cell  must  be  small  enough 
to  accommodate  steep  gradients  in  the  flow  properties  such  as  occur 
in  a  shock  wave.  Similarly  Atm  must  be  small  enough  to  avoid  gross 
changes  in  the  distribution  function  caused  by  either  the  collision 
or  the  convection  increment  over  the  time  interval  Atm.  There  is 
also  a  rough  lower  bound  on  the  cell  size  and  on  Atm.  This  is  because 
enough  collisions  must  be  computed  in  each  cell  during  Atm  to  assure 
that  the  computed  time  increment  is  fairly  close  to  Al^.  If  cell  k 
contains  n^  molecules  then  approximately  1/2  n^  collisions  will  be 
computed  per  collision  time.  Thus,  if  Atm  =  1/3  collision  times  and 
nk  =  30,  5  collisions  will  be  .computed  during  each  time  interval. 

The  applications  of  this  method  thusi  far  carried  out  have  used  30 
or  more  molecules  per  cell,  although  Bird  reports  in  one  of  his  papers 
[  6  ]  that  satisfactory  results  are  obtained  with  as  few  as  six  molecules 
per  cell. 


^°3.3  Storage  Requirements 

The  fast  access  storage  capacity  of  the  computer  limits 
the  number  of  molecules  that  can  be  used  in  the  simulation.  The 
programs  designed  by  the  author  have  used  500-700  Fortran  instructions 
leaving  space  for  close  to  1+000  5-coordinate  molecules  when  run  on 
the  University  of  Toronto  IBM  7091+  II  computer.  A  typical  program 
divides  naturally  into  an  initialization,  iteration,  and  output  section 
of  which  only  one  needs  to  be  in  fast  access  storage  at  one  time, 
thereby  reducing  further  the  storage  requirements  for  program. 

A  general  principle  is  that  the  larger  the  number  if  simulator 
molecules,  the  more  reliable  the  results.  In  the  case  of  time  independent 
problems  there  is  also  a  disadvantage  in  using  a  large  number  of  molecules 
because  computing  time  per  collision  time  increases  linearly  with  the 
number  of  molecules  and  a  fixed  number  of  collision  times  is  required 
to  achieve  convergence  to  the  steady  distribution.  This  unproductive 
computing  time  can  be  minimized  in  a  series  of  runs  where  the  parameters 
are  incremented  only  a  small  amount  each  run  by  using  the  final  sample 
from  one  run  as  the  initial  sample  for  the  next  run,  thereby  reducing 
the  convergence  time. 
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H.3.4  Computing  Time  Requirements 


For  purposes  of  estimating  the  actual  running  time  of 
a  given  program,  the  basic  parameters  are  tc,  the .time  required- 
to  compute  one  collision,  the  time  -required  to  move  one 
molecule  and  check  boundaries,  and  y  the  number  of  Ai^-s . per 
collision  time.  Since  approximately  1/2  H  collisions  must  be 
computed  per  collision  . time  (N  is  the  total  number  of -molecules 
and -each  collision  computation  accounts  for  two  molecules),  it  is 
easy  to  see  that  an  approximate  formula  for  the  running  time  of 
the  algorithm  is 

time  =  1/2  N  tc  +  H  n  =  1/2  N  (xc  +  2  p  Tip) 

Since  the  number  of  observations  of  a  given  variable  varies  as  the 
product  of  the  number  of  molecules  in  the  flowfield  and  the  number 
of  collision  times  involved,  it  follows  that  the  computing  time  per 
observation  is  constant  for  a  fixed  variable. 


Unfortunately,  the  best  data  available  in  the  literature 
gives  only  the  computing  time  necessary  to  compute  a  certain  number 
of  collisions.  Consequently  it  is  only  possible  to  determine  an 
approximate  value  for  the  combination  rc  +  2  y  tt  rather  than  for 
the  individual  parameters.  The  table  below  yields  a  rough  idea  of 
the  times  involved. 


Reference  foo. 


Computer  Used 


Computing  time/1000  collisions 


(3) 

(8) 

(12) 

Author 


Silliac,  U.  Sidney 
IBM  7040/32K 
cdc  66oo 
IBM  360/65 


60  x  1000/30000 
900  x  1000/2  x  106 
20  x  1000/6  x  105 
3  x  1000/10000 


=  2  min. 
=  .1*5  min. 

=  . 03  min . 

=  .3  min. 


U.3.5 


'he  Time ‘ Increment  Formula 


It  is  essential  in  a  simulator  of  the  type  we  have  been 
discussing  that  collisions  occur  at  the  right  rate  relative  to  the 
motion  of  the  molecules ,  and  that  each  pair  of  molecules  have  the 
correct  probability  of  undergoing  collision.  To  be  more  specific, 
oonsider  a  cell  %  in  position  space  of  volume  vol  (%)  and 
containing  nj,.  molecules.  This  means  that  contains  1/2  n^Cn^-l) 
distinct  pairs  of  molecules.  Let  g^j  denote  the  relative  speed  of 
molecule  i  and  molecule  j.  If  the  collision  cross-section  is  A, 
the  pair  must  be  found  in  a  volume  g. .AAtm  for  a  collision  to 
occur  during  time  At^.  Since  it  is  already  known  that  the  pair 
is  found  in  vol  (%),  the  probability  that  the  pair  i,j  undergoes 
collision  during  time  Atm  is 

g. .  A  At 

p  =  ^  J 

ij  voiTfy 
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and  the  probability  of  a  collision  is 


I 

i<J 


p«, 


where  tvere  are  L  =  1/2  %(%-!)  terms  in  the  sum.  Let  the  probabilities 
above  be  renumbered  from  1  to  L  so  that  given  i,j  there  exists  an  Z  such 
that  =  Pij ,  and  let  q£  =  l-p£.  This  means  that  the  probability  of 
exactly  m  collisions  is  the  sum  of  all  products  of -the  form 


**1  r2  ...  rL, 

where  r$,  =  pj,  or  q&  and  there  are  exactly  m  p&ts  in  .each  term. 

On  the  basis  of  the  foregoing  theory,  the  choice  of  collision 
partners  could  proceed  in  one  of  two  ways : 

A.  Accept  pair  Z  for  collision  with  probability  p^  and  repeat 
for  all  possible  pairs  Z  =  1,...,L.  Since  in  practice  there  are  about 
100  molecules  per  cell  yielding  roughly  50,000  pairs  of  which  only  a 
small  fraction  collide  during  time  Atm,  this  procedure .would  turn  out 
to  be  prohibitively  time  consuming. 

B.  Decide  the  number  of  collisions  at  the  start  by  choosing 

hjjj  and  then  choose  m 


This  method,  too,  is  unfeasible  because  it  requires  the  computation  of 
both  p^  and  hm  which  is  out  of  the  question. 

The  method  proposed  by  Bird  is  to  choose  Z  at  random  in  the 
range  (o,L)  and  r  from  the  uniform  distribution  on  (o,pta)  where  pba  is 
any  number  larger  than  all  the  p^-s.  The  integer  Z  is  accepted  if 
r  <  p;  otherwise  another  choice  is  made.  Once  Z  is  chosen,  a  counter 
c  is  advanced  by  1/Lp  and  the  collision  computation  is  complete  as  soon 
as  c  surpasses  1. 

It  will  be  shown  presently  that  this  procedure  gives  the 
correct  expected  number  of  collisions,  but  in  general  the  wrong 
distribution.  To  derive  the  expectation,  let  Z  =  1,,..,L  be 
random  variables  defined  by 


a  positive  integer  m  from  the  distribution 
successive  pairs  from  the  distribution 

L 


VM 


X 


z 


1  when  pair  Z  collides 
0  otherwise. 


Then  L 

s  =  I  X 
1  * 

gives  the  total  number  of  collisions  in  cell  k  during  time  At  , 


Furthermore 


L  L 

E(S)  =  Z  E(X  )  =  E  p 
1  £  1  & 

which  is  the  actual  expectation.  In  order  to  find,  the .Bird  expectation, 
observe  that  the  criterion  for  ending  the  collision  computation  is  a 
first  passage  in  a  random  walk  with  nonzero  mean  (see  Ref .2?2.10)jn 
such  a  process  the  expectation  of  the  time  to  first  passage  is  given 

fey 


E(t  at  time  of  first  passage) 
E(each  increment) 


Since  in  practice  the  increments  are  much  smaller  than  1, 
to  assume  that  the  numerator  will  be  very  close  to  1.  In 
compute  the  denominator,  note  that  the  Bird  procedure  x'o-** 
Z  yields  an  integer  from  the  distribution  (p„).  Then 

M 

E(increment)  =  E  p.  = 


&=1 


LZPo  2P0 


it  is  safe 
order  to 
choosing 


from  which  it  follows  that  on  the  average 

L 


choices  are  required  to  advance  the  parameter  c  by  one. unit.  In  order 
to  show  that  in  general  the  Bird  technique  does  not  produce  the  correct 
distribution  for  the  number  of  collisions  it  is  sufficient  to  display 
a  counter  example.  Consider  a  hypothetical  extreme  case  in  which 
P£  =  p  £  =  1,...,L.  Then  the  probability  that  m  collisions  occur  is 


h 

m 


m  L-m 
P  1  > 


(  )  being  the  binomial  coefficient,  whereas  the  Bird  method  would 
advance  c  by  1/Lp  for  each  collision  leading  to  the  probability 


P(m)= 


rl  if  m  =  [Lp] 
V.0  otherwise 


+1 


where  [Lp]  denotes  the  greatest  integer  less  than  or  equal  to  Lp. 


It. 3. 6  Boundary  Strategies 

There  are  two  commonly  occurring  boundary  situations  which 
must  be  handled  in  a  Monte  Carlo  simulation. 

a)  An  imaginary  boundary  in  a  drifting  Maxwellian  flow.  That 

is,  a  flow  whose  distribution  function  at  each  point  x  is  given  by 

_  _  2  3  2 

f(v)  =  n^(c  )~3  exp(-c  Ei=1  (v^-uj  )  (l) 
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where  the  vector  u  with  components  u^  i  =  1,2,3  specifies  the  velocity- 
of  the  flow,  c  is  a  constant  related  to  the  temperature  of  the  gas, 
and  nx  is  the  .local  number  density.  Any  molecules  reaching  such  a 
boundary  from  the  interior  of  the  region  are  simply  abandoned.  At 
the  end  of  each  time  interval  Atm,  it  is  necessary  to  introduce  a. 
number  of  molecules  representing,  those  that  have  crossed  the  boundary 
into  the  region  during  this  time .  To  be  specific ,  consider  -an  area 
element  6A  with  interior  normal  I  in  such  a  flow.  Let  (l,J,K)  be  an 
orthonormal  frame  such  that  J  lies  in  the  plane  of  I  and  the  velocity 
vector  u  of  the  flow.  Then  u  can  be  resolved  as  follows 


u  = 


"I1  +  V- 


(£) 


The  molecules  crossing  6A  into  the  .region  have  velocity  compr  '.ents  in 
the  I,J,K  directions  which  are  independently  distributed  and  have 
density  functions 

(l/c  it  1//2)  v.j.  expC-^Vj-UjJVc2) 

(l/c  7r  1/2)  exp(-(Vj-Uj)2/c2),  (3) 

and  (l/c  it  ^2)  exp(-(v^-u^)2/c2) 

respectively,  modulo  normalizing  factors.  The  total  flux  across  8 A 
is  given  by 

00 

<j)  =  n^  f  (v/c  it  1'2)  exp(-(v-uI)2/c2)dv  (1+) 

o 


In  nractice,  n  can  be  determined  from  the  formula 
x 


n  =  (XA  H,  (5) 

x  m 

-1/2 

where  2  X  is  the  local  mean  free  path,  and  Aj,,  is  the  cross-section 
per  simulator  molecule.  Then 


4>6AAtm 

gives  the  total  number  of  simulator  molecules  entering  the  region 
through  8k  during  time  Atm.  In  practice,  a  time  of  entry  is  chosen 
from  the  uniform  distribution  on  (0,Atm),  a  point  of  entry  is  chosen 
from  the  uniform  distribution  on  8k,  and  velocity  components  Vj,  Vj, 

Vj£  are  chosen  from  the  distribution  described  in  (3)  (Sec.  2.9.5). 

The  molecule  is  then  stored  with  the  position  coordinates  it  would 
have  at  the  end  of  the  time  interval  Atm. 

b)  Surface  interactions.  The  most  common  surface  interaction 

models  to  date  have  been  specular  and  diffuse  reflection.  In  the  case 
of  specular  reflection,  a  molecule  which  collides  with  a  boundary  during 
time  Atm,  has  its  velocity  vector  replaced  with  the  reflected  vector  at 
time  of  collision  and  moves  with  the  new  velocity  for  the  remainder  of 
the  time  interval.  The  reflected  velocity  v'  is  given  in  terms  of  the 
incident  velocity  v  by 


v"  =  v-2(r  •  1)1 
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where  I  is  the  normal  at  point  of  impact.  In  the  case  of  diffuse 
reflection;,  the  velocity  vector  is  replaced  by  one,  chosen  from  the 
distribution  (3)  (2.9.6)  in  which  u  ~-  0  and.  c  is  characteristic  of 
the  surface  temperature. 


4.3*7  Standard  Deviation  and  Limitations 

The  basic  theorem  relevant  to  the  discussion  of  the  standard 
deviation  irt  any  Monte  Carlo  method  asserts  that  the  standard  deviation 
of  an  average  of  independent  random  variables  decreases  as  the  inverse 
square  root  of  the  number  of  variables.  The  hypothesis  of  independence 
is  essential  for  it  is  even  intuitively  clear  that  in  the  .simulation  of 
a  sample  of  1000  molecules  for  one, collision  time  with  regular  observ¬ 
ations  of  the  distribution  function  during  this  time  there  comes  a  point 
after  which  nothing  more  is  gained  by  increasing  the  frequency  of 
observations.  This  is  because  eventually  successive  observations  will 
cease  to  be  independent  of  one  another.  This  issue  only  arises  of  course 
in  sampling  moments  of  the  distribution  function  at  regularly  spaced 
intervals.  It  does  not  arise  in  the  accumulation  of  information  on 
boundary  crossings  for  the  purposes  of  estimating  fluxes  because 
successive  boundary  crossings  are  automatically  independent. 

There  are  as  yet  no  studies  in  the  literature  of  how  the 
standard  deviation  of  a  typical  moment  depends  on  the  sampling  inter¬ 
val.  The  only  work  remotely  relevant  is  an  early  publication  by  Bird 
(3)  on  the  time  taken  for  the  distribution  function  to  relax  to  a 
Maxwellian  form.  It  was  of  the  order  of  3  collision  times.  This  is 
longer  than  the  sampling  interval  used  by  most  authors  in  the  field 
where  1/4  to  1  collision  times  has  appeared  to  be  satisfactory  for  most 
moments.  The  price  paid  for  the  wrong  sampling  interval  is  the  wrong 
theoretical  estimate  for  the  standard  deviations  and  wasted  computer 
time.  If  the  interval  is  too  short,  time  is  wasted  in  taking  needless 
samples,  and  if  too  long,  time  is  wasted  in  computation  between  samples. 

At  present  state  of  development.  Bird's  Monte  Carlo  method 
cannot  cheaply  yield  a  satisfactory  approximation  to  the  complete 
distribution  function  of  even  a  simple  practical  problem.  To  see  this, 
consider  the  flow  around  a  cylinder  at  Knudsen  number  1  and  make  the 
modest  demand  that  the  distribution  function  be  accumulated  on  a  grid 
having  10  points  in  each  direction,  leading  to  105  grid  points - 
altogether.  In  order  to  achieve  a  relative  standard  deviation  of  1/10 
at  each  grid  point,  roughly  100  observations  per  cell  are  required, 
which  is  107  observations  in  all.  On  the  assumption  of  a  computer 
memory  capacity  of  5  x  103  molecules,  107/5  x.  103  =  2  x  103  independent 
observations  on  the  flowfield  are  necessary.  On  the  basis  of  sampling 
interval  of  one  collision  time,  it  is  found  that  the  program  would  have 
to  run  for  2000  collision  times.  In  the  author's  experience,  the 
calculation  requires  about  one  minute  of  IBM  7094  II  or  IBM  360/65 
computer  time  per  collision  time.  Therefore  more  than  30  hours  of 
computer  time  would  be  needed. 

The  method  is  at  present  adequate  for  finding  certain  moments 
of  the  distribution  function.  In  order  to  obtain  a  density  field  with 
10$  relative  standard  deviation  at  102  grid  points  using  the  hypothetical 
program  of  the  last  paragraph  would  require  about  20  minutes  of  computer 
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time  which  is  well  within  reason.  Of  course,  a  relative  standard 
deviation  of  10$  is  rather  large,  leading  to  numbers  in  each  cell 
which  are  within  20%  of  their  mean  value  with  probability  0.95,  or 
within  30%  with  probability  0.997.  Results  of  such  coarseness  are 
usually  of  little  interest.  In  order  to  improve  them  by  a  factor  of 
■JW  '  requires  a  tenfold  increase  in  computing  effort,  producing 
a  relative  standard  deviation  of  0.036.  Then  one  is  assured  numbers 
within  7.2 %  of  the  mean  with  probability  0.95  or  9-5$  with  probability 
0.997-  Such  accuracy  is  acceptable  for  many  practical  applications. 


b.h  Variants  and  Refinements 


4,h.l  Weighted  Cells 

In  the  simulation  technique  described  thus  far  it  has  been 
understood  that  the  cross-section  per  molecule  was  the  same  for  each 
molecule  in  the  flowfield.  Another  way  of  looking  at  it  is  that  each 
simulator  molecule  represented  a  fixed  number  of  physical  molecules. 
There  are  many  situations  where  such  a  representation  is  unsatisfactory 
because  it  is  desirable  to  obtain  more  detail  about  the  flow  in  some 
subregion  cf  the  flowfield,  or  in  the  case  of  spherical  and  cylindrical 
coordinates,  to  obtain  information  near  the  origin  or  near  the  axis 
where  there  are  few  molecules  because  of  the  small  cell  volumes. 
Consequently  a  weighting  scheme  has  been  tried  by  both  G.  A.  Bird 
and  the  author. 

In  this  scheme  the  collision  computation  and  the  time  advance 
formula  remain  unchanged,  escept  that  the  cross-section  per  molecule 
becomes  dependent  on  the  cell.  Difficulties  arise,  however,  when  a 
molecule  crosses  a  cell  boundary  during  the  convection  phase  of  the 
iteration.  To  be  specific,  let  the  cross-section  per  molecule  in  cell 
k  be  Ajj  =  A/wj,.  k  =  1,,..,  K  where  w^  is  a  given  positive  cell  weight, 
and  A  is  a  constant  dependent  on  the  desired  total  number  of  simulator 
molecules  and  the  average  cross-section  density.  If  a  molecule  moves 
from  one  cell  to  another  having  twice  the  Vt-'ght,  it  is  reasonable  to 
expect  that  a  satisfactory  strategy  would  be  to  replace  it  with  two 
copies  of  itself.  The  dilemma  arises  when  the  molecule  moves  in  the 
opposite  direction.  A  resolution  which  has  been  found  satisfactory  by 
the  author  makes  use  of  a  random  integer  function  defined  as  follows: 

Intr(x)  =  [x]  +  i, 

where  [x]  is  the  greatest  integer  which  is  less  than  or  equal  to  x, 
and  i  is  a  random  variable  with  the  distribution 

j'O  with  probability  l-(x)-[x]) 

l  =  < 

U  with  probability  x~[x] 

Then,  when  a  molecule  moves  from  a  cell  with  weight  wa  to  a  cell  with 
weight  wg  ,  it  is  replaced  with 


intr(w0/vg) 

molecules  having  the  same  coordinates . 
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A  value  for  the  constant  A  compatible  with  the  molecule 
capacity  can  be. obtained  as  follows:  Suppose  it  is  given  that  the 
cross-section  density  ('l//5r  x  reciprocal  mean  free  path)  at  some  point 
and  time  in  the  flow  is  AQ.  A  useful  first  approximation  can  be  based 
on  the  assumption  that  the  .average  cross-section  density  in  the  whole 
flowfield  will  also  be  A D.  Then  cell  k  will  contain 

V°VAv‘ 

simulator  molecules,  where  vol^  is  the  volume  of  cell  k.  Then  the 
total  number  of  molecules  in  the  flowfield  will  be 


N 


av 


K 

SAw,  vol,  /A 


(&) 


where  N  depends  on  the  memory  capacity  of  the  computer.  Solving 
(6)  it  H  found  that 


K 

A  =  (A/N  )  S  w,  vol, 
av  ,  ,  k  Tc 
k=l 


4. If. 2  Mixtures  of  Gases 

The  Monte. Carlo  simulation  method  has  also  been  extended  by 
Bird  to  deal  with  binary  mixtures  of  gases.  Aside  from  minor  complic¬ 
ations  arising  from  the  3.abelling  of  the  molecules  in  computer  memory, 
the  biggest  modification  that  must  be  .made  is  in  the  time  advance 
procedure  in  order  to  ensure  that  the  correct  number  of  each  of  the 
four  possible  types  of  collision  is  computed.  Let  the  molecule  types 
be  numbered  1  and  2,  and  denote  the  collision  types  by  (l,l),  (l,2), 
(2,l),  and  (2,2).  The  basic  idea  is  to  compute  (l,l)  collisions  until 
the  time  parameter  has  been  advanced  by  1/4  At  ,  (1,2)  collisions  for 
another  1/4  At  etc . ,  using  the  same  time  increment  formula  as  before 
with  the  appropriate  cross-section  for  each  collision  type.  The 
somewhat  more  refined  procedure  adopted  by  Bird  entails  keeping  a 
separate  time  counter  for  each  collision  type  and  using  the  formula 


2  vol 


Aae\n6 


s 


to  increment  the  (a, 8)  counter.  At  each  step  in  the  collision  comput¬ 
ation,  a  collision  corresponding  to  the  counter  with  the  lowest  value 
is  computed  and  computation  ceases  when  each  counter  has  surpassed 
Atm,  The  formulas  for  computing  post  collision  velocities  for  molecules 
of  unequal  mass  are  commonly  found  in  textbooks  and  present  no  more  of 
a  problem  than  do  those  for  like  molecules. 
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It.  7  Definitions 


die” 


E(X) 

f  (ic,vst) 


(x*y) 


Volume  differential  in  3-dimensional  space 
(=  dx1dx2dx3,  where  xj,  x2,  X3  are  components 
of  the  vector  ,x) 

Expectation  of -the  random  variable  X. 

A  function  such  that -f(x,v,t)  dx::dV”  is  the 
average  number  -of  molecules  in  a  cell  of  6- 
dimensional  space  containing  the  point  .and  having 
6-dimensional  volume  "dx  <'  dv. 

Inner  product  of  the  vectors  x  and  y. 


(a, 3) 

& 

m 

[x] 


’Belongs  to’ 

v 

Interval  with  endpoints  a  and  3. 

-  P(P-l)  •••  ( L-m+1 ) 

1.2.'  ...  'm 


v 


the  greatest  integer  less  than  or  -equal  to  x. 
Length  of  the  vector  v. 


5.  MONTE  CARLO  EVALUATION  OF  THE  COLLISION  INTEGRAL 


5.1  Introduction 


The  last  method  to  he  discussed  in  this  review  has  been  developed 
almost  exclusively  by  Nordsieck,  Hicks,  and  their  coworkers  at 'the  Univ¬ 
ersity  of  Illinois,  and  described  reasonably  completely  in  a  series  of 
reports  of  the  Coordinated  Science  Laboratory.  Its  essential  features 
include  Nordsieck* s  proposal  from  as  early  as  1957  (2,  p.l)  for  choosing 
a  random  sample  of  points  in  the  region  over .which  the  collision  integral 
is  to  be  taken,  a  conventional  finite  difference  scheme  for  obtaining  a 
distribution  function  satisfying  the .Boltzman  equation  once  the  collision 
integral  is  known,  and  a  series  of  refinements  designed  to  overcome  approx¬ 
imation  and  random  errors. 

The  method  has  so  far  only  been  applied  to  problems  described 
by  the  hard  sphere  model  and  requiring  only  one  space  dimension  such  as 
the  .plane  normal  shock  wave  and  certain  relaxation  problems .  For  such 
problems  the  method  takes  about  one  hour  of  CDC  160U  time  for  the  iteration 
to  converge.  Each  iteration  yields  a  statistical  approximation  to  the 
iterate  one  would  obtain  with  an  exact  evaluation  of  the  collision  integral. 
The  approximation  has. a  standard  deviation  of  about  3 %  at  each  of  10  x  226 
=  2260  histogram  bins  (10  position  cells  and  226  velocity  bins). 

Extension  of  the  method  to  problems  requiring  higher  space 
dimension  seems  difficult  because  of  computer  storage  and  time  require¬ 
ments  and  also  because  of  the  problems  involved  in  devising  a  successful 
finite  difference  scheme  in  more  than  one  dimension. 

The  discussion  presented  here  will  also  deal  exclusively  with 
problems  requiring  one  space  dimension  and  will  treat  as  an  example 
Hick's  and  his  coworkers'  calculations  of  the  plane  normal  shock  wave. 


5.2  Theoretical  Preliminaries 


5.2.1  Nordsieck 's  Units 


Hicks  and  his  coworkers  have  made  use  of  special  units  .derived 
basically  from  the  standard  procedure  for  nondimensionalizing  the  Boltzman 
equation.  The  units  are  based  on  conditions  in  some  reference  region  of 
the  flow  (in  the  shock  wave  problem,  the  far  upstream  conditions  are 
convenient) . 

If  bm,  m  are  the  molecular  diameter  and  mass  respectively,  k 
is  Boltzman 's  constant,  and  T,n  are  the  reference  temperature  and  number 
density  respectively,  the  Nordsieck  units  for  length,  velocity  and  time 
are  respectively  Z  =  l/2trn  bm2,  c  =  (2  kT/m )1'2,  and  Z/c,  and  the 
distribution  function  is  given  in  terms  of  n/c3. 

If  the  substitutions 

v=Ax,  v  =cv  ,  v  "=cv  ",  v^cvl",  t=U/c)t,  f=(n/c3)f , 

X  A  X  X 

f*=(n/c3)f",  F=(n/c3)F,  and  F"=(n/c3)F" 
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»v 

are  made  in  BoltzmanA-s  'equati'on*  tlientdfter  a3h6ppihg.'the  -  -s^/ it (becomes 
(See- (di);rl<2. .5 't '' 


In  their  calculations l'  Nordsieck,  Hicks,  ef  al  make  use  of  a  further 
velocity  scale  change  in  order  to  allow  the  use  of  velocity  bins  of  a 
size  convenient  for  computation.  Since  this  change  is  not  necessary  for 
explaining  the  essential  features  of  the  method,  it  will  he  omitted  here. 


5.2.2  The  Discrete  Analogue 


In  the  computational  algorithm,  both  velocity  and  position 
space  are  quantized.  Furthermore,  in  the  case  of  velocity  space,  some 
means  is  needed  to  represent  infinite  two  dimensional  space  using  only 
a  finite  number  of  cells.  This  last  problem  is  handled  in  the  Nordsieck 
-Hicks  method  by  restricting  attention  to  a  semi-circular  region  in  (v  , 
Vl )  space  which  is  still  large  enough  to  include  99%  of  the  molecules 
under  conditions  expected  in  the  flow. 


With  regard  to  quantization,  Nordsieck  and  Hicks  have  reported 
good  results  with  a  non-uniform  net  of  about  10  points  in  the  x-direction 
and  an  array  of  226  square -bins  in  (vx,  v^ ) -space  covering  the  semicircular 
region  with  2k  bins  along  the  diameter.  The  net  in  the  x-direction  is, 
non-uniform  because  if  varying  density  gradients.  In  fact,  in  their 
later  work  on  the  shock  problem,  Hicks  and  Smith  found  that  x  was  not  a 
good  variable  to  work  with  because  of  the  problems  in  making  the  .x-cells 
fine  enough  to  take  care  of  the  large  flow  gradients  in  the  shock  itself. 
After  some  experimentation  with  various  monotonic  transforms  of  x,  it 
was  found  that  the  best  variable  to  use  was  the  density  n  itself. 


Since  it  is  beyond  the  scope  of  the  present  work  to  go  into  the 
details  of  a  particular  solution  it  will  henceforth  be  assumed  for  sake 
of  simplicity  that  there  are  position  cells  with  partition  points  xj 
j=l,  ...  ,  jm,  and  s  square  velocity  bins  with  midpoints '(vxs»  vxS), 
s=l,  ...  ,  sm.  The  values  of  the  distribution  function  will  be  corres¬ 
pondingly  indexed  as  f(,j,s).  Then,  if 

S(s,s^,n)  and  S'lsjS^n) 


denote  the  cells  containing  the  post-collision  velocities  arising  out 
of  a  collision  between  molecules  having  velocities  v$  and  vs^  with 
collision  vector  n,  the  discrete  approximations  to  F  and  F'  are 

f(j,S)  and  f(j,S') 


respectively. 

The  discrete  approximation  to  (l)  suggested  by  Nordsieck  for 
the  time  independent  problem  where 
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is  then 


Vl"x«5 


=  1/2  [(A-Bf)jjS+(A-Bf)J+1^l 


where  (a-bf ) j  jS  is  an  approximation  to.  the  collision  integral  given  by 
s  q_  r  , 

( „  —  r  ®  a-.  a  r  a  ^  v  ^  _  „  I  r  .  (Tr  «  W  Y 


(a-bf)  =  l  m  v  ^Av  ..Av  ..  l  “  Ato'  E  m  —  |n  *(v  -v  „)|X 
■>’s  s'=l  is  xs  is  4=1  fl  r=lrm  r  s  s 


X  [f(j,S)f(j,S')-f(3,s)f(3,s')l 

where  n  r=l,  .  .  .•  *  r  is  an  equidistant  array  of  points  on  the  unit 
sphere  and  m 

qL 

T,  A  ip'  =2ir. 

q=l  f<L 


The  evaluation  of  the  discrete  approximation  to  the  collision 
integral  in  (6)  (a  sum  repeated  jmsm  times)  is  still  beyond  the  capabilities 
of  the  present-day  computer.  To  date,  the  only  resolution  of  the  diffic¬ 
ulty  is  the  Monte  Carlo  method  proposed  by  Nordsieck.  Nordsieck's  method 
is  based  on  conventional  Monte  Carlo  evaluation  of  the  sum  approximating 
the  collision  integral  and  not -of  the  , integral  itself.  The  essence  of  the 
procedure  is  as  follows:  In  order  to  evaluate  the  sum  at  fixed  (j,s),  a 
random  sample  of  N  values  of  the  quadruple 

<v>  V,  ?'■  3> 

is  chosen,  where  (v^,  vj/)  is  uniformly  distributed  over  the  cell  centres 
of  velocity  space,  "  is  uniformly  distributed  over  an  equally  spaced  set 
of  points  in  (0,  2m),  and  n  is  uniformly  distributed  over  an  equally 
spaced  set  of  points  on  the  surface  of  the  unit  sphere.  The  Monte  tarlo 
approximation  to  (a-bf)jjS,is  then  given  by 

(3) 

The  method  as  it  is  presently  used  differs  from  the  above  in 
a  way  which  takes  advantage  of  certain  symmetry, ' exchange  and  reflection 
properties  of  the  collision  process  in  order  to  expand  the  effective  size 
of  the  sample  at  a  considerable  saving  in  work.  In  order  to  make  this 
possible  it  is  necessary  to  evaluate  times  the  collision  integral 

/“  r00  r ^  a~ 

dv'J  VjVj/dVj/J  g-  |n-(v  -v)l(FF'-ff')  (4) 
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instead  of  the  collision  integral  itself. 


The  details  are  as  follows:  two  three-dimensional  precollision 
vectors  v,  v'  are  chosen  instead  of  only  one.  Each  choice  entails  two 
steps . 

A.  (vx,v,_)  is  chosen  from  among  the  cell  centres  in  velocity  space. 

B.  •  (p  is  chosen  from  an  equally  spaced  set  of  points  in  (0,27r). 


The  unit  vector  n  is  then  chosen  at  random  from  an  equally  spaced  set 
of  points -on  the  surface  of  the  unit  sphere ,  and  the  post-collision 
velocities  V,  V'  are  computed.  It  is  convenient  to  call  a  sample  of 
3-vector  5-tuples 


admi s s ible '( AS )  i f  the  individual  vectors  have  the  properties  just 
Then  the  -expression 


described. 


226 

N 


Ss,  fe yLV|5- (5-5')| 


(5) 


where  s,  s',  S,  S'  are  cell  numbers  of  v,  v',  V,  V'  respectively  is  an 
approximation  to , ( 4) . 


The  advantages  of  this  approach  become  apparent  when  it  is 
observed  that  if 


((v,  v',  n »  V,  V')j 


is  admissible  then  so  is 


[(v%  v,  n,.  V',  v| 


(6) 


Furthermore,  if  the  velocity  space  is  symmetrical  with  respect  to 
reflection  about  v  =0  and  the  subscript  R  denotes  this  reflection,  then 

X 


{^R5  V  R*  n>  V  ^  r| 


(7) 


is  also  admissible. 


Finally,  if  v,  v',  n  V,  V'. denotes  a  collision,  then 
V,  V',  -n  ->  v,  v'  also  denotes  a  collision  called  the  inverse  collision. 
In  addition,  since  the  Jacobian  of  the  transformation 


(v,  v')  -v  (V,  ¥') 


has  absolute  value  1  for.  each  n,  it  follows  that  if  each  of  v,  v'  are 
uniformly  distributed  over  some  region  of  velocity  space,  V,  V'  will  be 
uniformly  distributed  over  the  image  region.  The  image  region -in  Nordsieck’s 
case  is  slightly  different  from  the  original  region  so  that  collisions 
yielding  V,  V'  outside  if  the  original  region  must  be  discarded.  The 
discarded  fraction  is  quite  small  (less  than  0.l6  [9])  and  occurs  on  the 
fringe  of  the  finite  velocity  space  where  the  distribution  function  f  is 
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(8) 


also  very  small.  It  is  expected  that  the  error  arising  from  this 
discrepancy  is  insignificant  so  that 


can  also  he  considered  admissible. 

In  view  of  the  preceding  discussion,  it  is  clear  that  starting 
any  admissible  sample,  seven  new  admissible  samples  can  be  constructed 
by  utilizing  (6),  (7),  and  (8),  and, the  composite  sample  will  be  eight 
times  as  large  as  the  original  one.  It  is  equally  clear,  however,  that 
the  composite  sample  will  not  consist  of  statistically  independent 
choices,,  although  it  will  be  a  fair  sample .for  purposes  of  evaluating 
the  collision  integral.  This  lack  of  independence  does  not  receive 
comment  in  the  literature  on  Nordsieck's  method,  but  it  seems  likely 
that  if  anything  it  is  beneficial  in  .reducing  the  variance  of  the 
estimate  for  the  collision  integral.  The  actual  strategy  for.  incorp¬ 
orating  the  composite  sample  into  the  collision  integral  estimate  is 
described  as  part  of  the  algorithm. 


|{V,  V%  -n,  v,  v')j 


5.2.4  The  Corrections 

Hicks,  Yen,  and  Reilly  have  published  a  report  dealing  exclusively 
with  various  techniques  of  error  reduction  and  giving  the  formulas  which 
they  have  found  through  experience  to  be  most  effective.  In  most  cases, 
theoretical  support  for  the  correction  procedures  is  lacking. 


Errors  enter  into  the  collision  integral  estimate  from  two 
sources:  the  replacement  of  the  collision  integral  by  a  discrete  approx¬ 
imation  and  from  the  statistical  fluctuations  inherent  in  the  Monte  Carlo 
method.  It  was  found  very  early  that  if  nothing  was  done  to  normalize  f 
after  each  iteration,  a  systematic  drift  in  the  values  of  f  would  lead  to 
an  eventual  violation  of  the  laws  of  conservation  of  mass,  momentum,  and 
energy.  This  difficulty  is  presently  handled  by  the  least  squares  (LS) 
correction. 


In  the  (LS)  correction,  corrected  values  ps  s=l,  ...,  sm,  of 
the  collision  integral  are  chosen  in  such  a  way  as  to  minimize  the  sum  of 
squares 


m 


I  [(A-Bf)  -  p  ]2/p2 

15  5  S 


subject  to  constraints  arising  from  the  conservation  laws.  In  the  shock 
problem,  mass  momentum,  and  energy  fluxes  are  conserved  leading  respect¬ 
ively  to  the  constraint  equations 


Z  p  v,  v 
■'s  J-s  xs 


Z  p  v.  v  “ 
is  xs 
s 


const 

const 


Z  p  v  v  (v  2+v,  2)=  const 
rSiS  xs  xs  is 
s 
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The  correction 


A  =  p  -  (A-Bf)  ,  s=l,  s,  (9) 

s  s  s  hi 

is  split  evenly  "between  As  and  (Bf)s  except  in  the  case  of  velocity  bins 
where  one  of  the  corrected  values  would  be  negative,  when  it  is  split 
in, such  a  way  as -to  leave  both  terms  positive. 

An  attempt  is  made  to  control  the  error  arising  from  statistical 
fluctuations  by  the  Maxwell-Bolt zman  (MB)  correction.  In  this  correction, 
the  values  of  the  collision  integral  are  calculated  in  such  a  way  that  the 
collision  integral  would  be  zero  for  a  gas  in  equilibrium  which  has  the 
same  values  of  density,  temperature,  and  velocity  as  exist  at  the  given 
point  in  the .non-equilibrium  gas.  In  carrying  out  this  .correction,  the 
required  moments  of  f  must  first  be  computed  in  order  to  obtain  for  each 
position  station  j  the  local  density  n j ,  the  mean  speed  uj  and  a  measure 
of  the  deviation  from  the  mean  c j .  Let  fMB  denote  the  Maxwell-Bolt zman 
distribution  function  based  on  these  parameters,  i.e. 

WVV  =  2nj7r"1/2Cj"3viexp-[vx-uJ)2  +  vx2]/c2 

The  correction  corresponding  to  a  collision  (v,  v%  n,  V,  V")  is  given 
by 

A(v,v;5)  =  1/2  vxvJLVn-(v"-v)|[fMB(V)fMB(V')-f(v)f(v)]  (10) 

Then  each  time  during  the  statistical  computation  of  the  collision 
integral  that  (v^A)(v)  and  (v^BfHv)  are  incremented -as  a  consequence  of 
such  a  collision,  each  is  immediately  corrected  according  to  the  formulas 

(V|A)c(v)  =  (VjA) (v)  -  A(v,v',w) 

(ll) 

(vLBf)c(v)  =  (vxBf)(v)  +  A(v,v',w) 


5.2.5  The  Finite  Difference  Scheme 

Nordsieck's  proposal  for  a  stable  finite  difference  scheme 
(one  in  which  old  errors  decrease  from  step  to  step)  makes  essential 
use  of  the  fact  that  the  loss  term  Bf  in  the  collision  integral  is 
proportional  to  f  and  that  Boltzman's  equation  in  the  one  dimensional 
time  independent  case  is  (see  1.3.1).  . 

|f  =  f-  (A-Bf)  (12) 

X 

The  elementary  theory  of  first  order  linear  ordinary  differential 
equations  says  that  a  perturbation  of  the  initial  value  for  (12)  will 
decay  as  x  increases  if  vx  >  t>. 

Consequently,  the  direction  of  integration  in  any  finite 
difference  scheme  for  this  equation  should  be  to  the  left  if  fx  <  0 

and  to  the  right  if  vx  >  0.  In  the  case  when  x.-x.  =  Sx  =  const, 

0  0—1 
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a  typical  recursion  formula  for  integration  to  the  right  based  upon  the 
approximation  to  the  derivative  given  in  (6)  is  given  by 


f  =  f 
j+1  j  1+B 


*  *  r 


3+1 


l+B 


6+1 


where 


A'  =  A.  6x/2v  ,  B'  =  B.  6x/2v  ,  j=l,...,j: 

6  J  x5  j  o  •  x’  m 


This  simple  formula  did  not  prove  satisfactory  for  shock  wave 
studies  because  of  large  gradients  in  flow  properties  near  the  .centre 
of  the  shock.  Recent  work  by  Hicks  et  al  on  this  problem  makes  use  of 
n  as  an  independent  variable  instead  of  x.  To  write  Bolt zman ' s  equation 
in  terms  of  n  requires  knowing  the  function  dn,  In  practice  it  is 

cbr' 

obtained  at  each  step  of  the  integration  procedure  by  using  a  discrete 
version  of  the  formula 


dn 

dx 


(13)  , 


where  (A-Bf )^c  is  the  Monte  Carlo  approximation  to  the  collision  integral. 


5.3  The  Algorithm 

5.3.1  Storage  Requirements  and  Initialization 

The  major  blocks  of  storage  in  this  method  must  be  assigned 
to  storing  the  current  iterate  of  the  distribution  function  f  and  the 
two  parts  A,  Bf,  of  the  Monte  Carlo  estimate  of  the  collision  integral. 

If  jm  stations  are  used  along  the  x-axis  and  sm  bins  are  used  in  velocity 
space,  then  jmsm  locations  are  required  for  each  of  these  functions.  The 
values  of  A  and  Bf  must  be  stored  separately  (instead  of  the  difference 
A-Bf)  because  they  are  required  in  the  stable  difference  scheme  used  in 
integrating  the  sm  ordinary  differential  equations  (one  for  each  velocity 
bin) . 


In  addition  to  these  major  blocks  of  storage,  general  workspace 
is  required  as  well  as  space  to  store  certain  function  values  which  are 
used  repeatedly  and  which  would  take  an  excessive  amount  of  time  to 
recalculate  each  time  they  are  used. 

The  main  part  of  initialization  entails  giving  initial  values 
to  the  distribution  function  f,  that  is  specifying  the  zeroth  iterate. 

Not  unexpectedly,  best  results  are  obtained  when  the  initial  approxim¬ 
ation  is  fairly  close  to  the  correct  steady  state  distribution.  In  the 
case  of  Hicks  and  coworkers  shock  wave  studies,  the  initial  distribution 
was  taken  to  be  Mott-Smith. 


5.3.2  Steps  in  One  Iteration 

A.  Monte  Carlo  estimate  for  the  collision  integral.  As  has  already 

been  mentioned  in  5-2.3,  it  is  convenient  to  first  estimate  v^A  and  Vj.Bf 
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and  then  to  divide  by  v^  and  vxf  respectively,  in  order  to  obtain- A 
and  B.  Also,  the  Maxwell-Bolt zman  correction,  if  it  is  used,  must  be 
incorporated  here.  Finally,  it  is  worth  pointing  out  that  the  estimation 
procedure  about  to  be  described  must  be  repeated  times,  once  for  each 
value  of  the  position  variable. 

Before  commencing  accumulation  of  the  integrals,  it  is  necessary 
to  zero  the  storage  for  vxA  and  vxBf,  and  to  compute  the  factor 

w  =  2ir(i  /N  )  x  (area  of  velocity  space), 
m  s 

where  Ns  is  the -number .of  collisions  in  the  collision  sample,  which 
according  to  the  experience  of  Hicks  must  be  about  213  or  larger.  Thus 
the  sampling  procedure  about  to  be  described  is  repeated  Ns  times. 

In  this  procedure  a  velocity  pair  (vx»  vx)  is  chosen  at  random 
from  the  discrete  two  dimensional  velocity  space,  and  an  angle  (p  is 
chosen  from  an  equidistant  set  of  values  in  (0,  2ir).  It  must  be  remarked  - 
that  since  the  actual  value  is  not  needed  in  the  computation,  it  is  time¬ 
saving  to  store  only  the  sets  of  values  sin  (p  and  cos  (p  and  to  choose 
the  appropriate  one  of  these  at  random  when  the  need  arises.  In  any  case, 
the  triple  (vx,  vx ,  to  )  defines  a  velocity  vector  henceforth  denoted  by 
v..  The  procedure  is  thus  repeated  to  choose  another  velocity  vector  v" 
independently  of  the  first.  Finally,  the  collision  vector  w  is  chosen 
at  random  from  an  equally  spaced  set  of  points  on  the  surface  of  the  unit 
sphere . 


It  is  now  possible  to  compute  the  function 
X  =  Vjv/lw-vj^v,  =  v'-rv) 

and  the  post-collision  velocities  V,  V"  (lV.l.iL,Eq.l,2) .  Let  s,s',S,S' 
denote  the  cell  numbers  of  the  velocities  v,v',V,V'  respectively.  Finally, 
if  v  is  any  vector,  let  vp  denote  the  reflection  of  v  in  the  plane  x=0,  and 
let  sp,s"r,Sr,S'r  denote  the  reflected  subscripts  corresponding  to  s,s",S,S'' 
respectively.  One  further  piece  of  notation  is  convenient  in  writing  down 
the  subsequent  formulas.  That  is: 

f  =  f(j,s),  f'  =  f(j,s') 

F  -  f ( j ,S) ,  F"  =  f(jiS'), 


with  fR, 

from  the 


fp"* ,  Fp,  Fr'  defined  correspondingly. 

Eight  contributions  to  the  collision  integral  are  then  computed 
triple  (v,v%n)  by  adding 

wXFF  to  Vj_Als),  v^A(s') ,  Vj_Bf  ( S ) ,  vxBf{S'), 

wxff  to  vxA(S),  vxA(S") ,  vxBf (s) ,  v,.Bf\s'), 

VXFRFR'  to  vxA(sr),  vxA(s.^),  vxBf(SR),  vxBf(SR'), 

wxfRfR'  to  vxA(Sr),  vxA(Sr'),  vxBf(sR),  vxBf',sR'). 


7? 


At  the  same  time  the  MB  correction  can  be  applied  according-  to  the  formulas 

(vlA)c(s)  =  VjA(s)  -  A(s,s%w);  (vLbf)c(s).  =  Vjbf(s)  +  A(s,s',w), 

(v,A)  (s')  =  v.A(s')  -  A(s',s,n);  etc. 
c 

As  soon  as  Ns  sets  of  eight  contributions  apiece  have  been 
computed,  the  whole  Monte  Carlo  collision  estimate  is  repeated  at  the  next 
position  station.  Hicks  and  his  coworkers  found  that  the  properties  of 
the  general  iteration  procedure  were  improved  if  the  statistical  fluc¬ 
tuations  between  position  stations  were  eliminated  by  making  use  of  the 
same  collision  sample  at  each  position  station.  This  can  be  simply 
achieved  by  reinitializing  the  random  number  generator  to  the  .same 
value  at  the  start  of  each  integral  computation.  When  all  the  collision 
integrals  have  been  computed,  the  program  proceeds  to  B. 

B.  The  least  squares  correction.  The  LS  correction  is  the  smallest 
modification  of  the  sm  values  of  the  collision  integral  consistent  with 
conservation  of  mass,  momentum,  and  energy.  The  correction  is  made  in 
accordance  with  the  formulas  given  in  2db, 

C.  The  numerical  integration  phase.  Once  the  corrections  have 
been  made,  values  of  a  and  b  suitable  for  use  in  the  numerical  integration 
scheme  can  be  calculated.  The  numerical  integration  itself  involves  a 
solution  of  Sm  uncoupled  ordinary  differential  equations,  one  for  each 
fibre  sm  =  const  in  the  discrete  phase  space.  The  relevant  recursion 
relations  are  discussed  in  2„l4o  The  jmSm  values  of  f  obtained  as  a 
result  of  the  integrations  constitute  the  next  iterate. 

D. _  Calculation  of  moments.  Some  or  all  parts  of  this  subprogram 

are  generally  bypassed  until  the  final  iteration.  Normally,  only  enough 
information  needs  to  be  output  from  the  intermediate  iterations  to  give 
evidence  that  the  iterations  are  in  fact  converging.  It  would  be  too 
lengthy  to  give  a  complete  list  of  what  type  of  information  about  the 
distribution  inunction  has  or  should  be  calculated.  Nordsieck,  Hicks,  et 
al  have  computed  moments  up  to  the  fourth  in  the  shock  problem,  graphical 
displays  of  the  i^stribution  function,  and  estimates  of  standard  deviation. 

E.  Exit  criterion.  The  common  criteria!  for  convergence  of  an 
iteration  scheme  are  applicable  here.  The  simplest,  of  course,  is  a 
test  that  subsequent  Iterations  are  close  enough  together.  It  has  been 
found  that,  starting  with  the  Mott-Smith  distribution  in  the  shock 
problem,  around  12  iterations  give  satisfactory  convergence. 

5.k  Comments 

The  only  aspect  of  this  method  which  is  rigorously  understood 
from  the  standpoint  of  theory  is  the  Monte  Carlo  estimate  of  the  collision 
integral.  Here,  standard  results  on  evaluating  integrals  by  means  of  a 
random  sample  of  points  m  the  domain  of  integration  apply.  The  fact 
that  the  eight  variants  arising  from  each  collision  are  dependent  has 
not  been  treated  theoretically  but  presumably  comes  under  the  heading 
of  the  theory  of  antithetic  variates.  In  addition,  m  Nordsieck1 s  method 
it  is  not  the  collision  integral  itself  that  is  estimated,  but  only  a 
discrete  approximation.  The  closeness  of  this  approximation  has  only  been 
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investigated  in  an  empirical  sort . of  way.  In  the  numerical  integration 
phase,  the 'recently  opened  study  of  differential  and  difference  equations 
with  coefficients  chosen  from  some  distribution  could  conceivably  throw 
light  on  possible  inaccuracies  in  Nordsieek's  solution  of  the  226  diff¬ 
erence  equations  whose. coefficients, have  been  obtained. from  the  collision 
•integral  estimate. 

Hicks  and  his  coworkers  did  initially  run  into  difficulties 
on  the  numerical  integration  phase  which  they  first  tried  to  surmount 
by. smoothing. .the  values  of  a  and  bf  after  they  were. computed,  and  most  . 
recently  by  using  a  fixed  collision  sample  for  each  position  station. 

The. latter  .has-succeeded  in  stabilizing  the  calculation,  but  does  not 
help  in  obtaining  a  measure  of  how  far  the  calculated  solution  is  to  the. 
•actual  one.  A  great  deal  of  work  remains  to  be  done. 

■  At  the  moment ,  faith  in  the  method  must  rest  on  the  exhaustive 
empirical-  error-  studies  that  have  been  done  for  each  of  the  problems 
undertaken  .by  the  Chicago  team,  and  described  in  detail  in  their  CSL  •• 

.  reports .  The  care  with  which  this  work  has  been  done  evokes  confidence 
in  the  actual  figures  that  have  been  published.  Nevertheless,  the 
sparseness  of  theory  dealing  with  the  method,  the  intricacy  of  the 
corrections  which  seem  essential  to  its  success,  and  the  time  required 
to  both  develop  the  program  for  a  particular  problem  and  for  the  comput¬ 
ation  itself,  limit  its  usefulness  as  a  general  purpose  tool  for  solving 
kinetic  theory  problems.  Finally,  the  method  has  been  developed  for 
zero  or  one  space  dimension,  whereas  many  of  the  more  interesting 
problems  require  at  least  two. 
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5.6  List  of  Symbols 
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First  .Occurrence 

a 

5.2.2 

A,As 

5. 2. U 

V  (A-Bf)MC 

5.2.5 

b 

m 

5.2.1 

b 

5.3.2 

E  (Bf) 

s 

5.2.4 

bj 

5.2.5 

c 

5.2.1 

cd 

5. 2. 4 

r,  f,  fj  i',  f,  f%  p,  F' 

5.2.1 

f(d,s),  r(3,sO 

5.2.2 

fMB 

5.2.4 

fj 

5.2.5 

fR’  fR5  FR?  ?R 

5.3.2 

J]  5 

5.2.2 

A 

5.2.1 

m , 

5.2.1 

n,  n 

5.2.1 

n 

r 

5.2.2 

N 

s 

5.3.2 

R 

5.2.3 

s,  s  ,  S(s,s',n),  S(s,  s',  n) 

5.2.2 

T 

5.2.1 

Uj 

5.2.4 

Roman- 


Roman-  ■ 
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v. ,  v  ,  v',  v',  V  .  vf 
x’  x?  x’  x*  * 

5.2.1 

v  >  v  ,  V  ,  V" 
xs’  s’  s’  s 

5.2.2 

(vAa)(v),  (vxtf)(v),  (vAa)  (v),  (v.bf)  (v) 

c  c 

5.2.U 

VR 

5,3.2 

A 

w  ■ 

5.3.2 

X,  X 

5.2.1 

XJ 

5.2.2 

Greek 
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Avxs'’  AvV  % 

5.2.2 

A(v.,  v',  n) 
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